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Automated Parameter Selection for Total Variation 
Minimization in Image Restoration 

Andreas Langer 


Abstract Algorithms for automatically selecting a scalar 
or locally varying regularization parameter for total 
variation models with an L'^-data fidelity term, r S 
{1,2}, are presented. The automated selection of the 
regularization parameter is based on the discrepancy 
principle, whereby in each iteration a total variation 
model has to be minimized. In the case of a locally vary¬ 
ing parameter this amounts to solve a multi-scale total 
variation minimization problem. For solving the consti¬ 
tuted multi-scale total variation model convergent first 
and second order methods are introduced and analyzed. 
Numerical experiments for image denoising and image 
deblurring show the efficiency, the competitiveness, and 
the performance of the proposed fully automated scalar 
and locally varying parameter selection algorithms. 

Keywords Total variation minimization • Locally 
dependent regularization parameter • Automated 
parameter selection • L^-data fidelity • L^-data fidelity • 
Discrepancy principle • Constrained/unconstrained 
problem • Gaussian noise • Impulse noise 

1 Introduction 

Observed images are often contaminated by noise and 
may be additionally distorted by some measurement 
device. Then the obtained data g can be described as 

g =N{Tu), 

where u is the unknown original image, T is a linear 
bounded operator modeling the image-formation de¬ 
vice, and J\f represents noise. In this paper, we consider 
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images which are contaminated either by white Gaus¬ 
sian noise or impulse noise. While for white Gaussian 
noise the degraded image g is obtained as 

g = Tu + r], 

where the noise rj is oscillatory with zero mean and stan¬ 
dard deviation cr, there are two main models for impulse 
noise, that are widely used in a variety of applications, 
namely salt-and-pepper noise and random-valued im¬ 
pulse noise. We assume that Tu is in the dynamic range 
[0,1], i.e., 0 < Tu < 1, then in the presence of salt-and- 
pepper noise the observation g is given by 

{ 0 with probability ri G [0,1), 

1 with probability r 2 G [0,1), (l.I) 

Tu{x) with probability 1 — ri — r 2 , 

with I — ri — r 2 > 0. If the image is contaminated by 
random-valued impulse noise, then g is described as 

p with probability rG[0,l), (12) 

Tu{x) with probability 1 — r, 

where p is a uniformly distributed random variable in 
the image intensity range [0,1]. 

The recovery of u from the given degraded image g 
is an ill-posed inverse problem and thus regularization 
techniques are required to restore the unknown image 
|41j . A good approximation of u may be obtained by 
solving a minimization problem of the type 

min'H{u] g) + alZ{u), (1-3) 

U 

where 'H{.\ g) represents a data hdelity term, which en¬ 
forces the consistency between the recovered and mea¬ 
sured image, TZ is an appropriate filter or regulariza¬ 
tion term, which prevents over-fitting, and a > 0 is 
a regularization parameter weighting the importance 
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of the two terms. We aim at reconstructions in which 
edges and discontinuities are preserved. For this pur¬ 
pose we use the total variation as a regularization term, 
first proposed in |81] for image denoising. Hence, here 
and in the remaining of the paper we choose 7 ?.(m) = 
\Du\, where \Du\ denotes the total variation of 
u in 17; see wm for more details. However, we note 
that other regularization terms, such as the total gen¬ 
eralized variation m, the non-local total variation |60] , 
the Mumford-Shah regularizer E], or higher order reg- 
ularizers (see e.g. [77] and references therein) might be 
used as well. 


noise models, as Poisson noise EH, multiplicative noise 
m, Rician noise gS). For images which are simultane¬ 
ously contaminated by Gaussian and impulse noise |15j 
a combined L^-L^-data fidelity term has been recently 
suggested and demonstrated to work satisfactory [^ . 
However, in this paper, we concentrate on images de¬ 
graded by only one type of noise, i.e., either Gaussian 
noise or one type of impulse noise, and perhaps addi¬ 
tionally corrupted by some measurement device. 


1.2 Choice of the scalar regularization parameter 


1.1 Choice of the fidelity term 


For the reconstruction of such images the proper choice 
of a in m and A in (US is delicate; cf. Fig. 11.11 In 
particular, large a and small A, which lead to an over¬ 
smoothed reconstruction, not only remove noise but 
also eliminate details in images. On the other hand, 
small a and large A lead to solutions which fit the 


The choice of 'H typically depends on the type of noise 
contamination. For images corrupted by Gaussian noise 
a quadratic L^-data fidelity term is typically chosen and 

has been successfully used; see for example fTOlfroilTni , , i i r . . . , 

i7TTirT7nrT7TiiT77TirT7TirTFirr>7irT7TirT7ni77^ T • 1 1-1 given data properly but therefore retain noise in ho- 

[7ill28U291l301l32U35ll47ll721l761l911l93 . In this approach, which^ tt 

r i rmr i i • - • mogcneous regions. Hence a good reconstruction can 

we refer to as the L^-1 V model, the image u is recov- , , . , , , . , . , ^ i , 

, „ , 1 1 i 1 1 ■ be obtained by choosing a and respectively A such that 

ered from the observed data g by solving . ° ^ i m 

a good compromise ot the aforementioned effects are 

made. There are several ways of how to select a in 
(HI) and equivalently A in dH, such as manually by 
the trial-and-error method, the unbiased predictive risk 
where denotes the space offunctions with bounded estimator method (UPRE) pilHTj . the Stein unbiased 

T i r\ t l—C 1/ f I / \ iF n 11 r i F j / / A / f J \ _ 


mm — 
u^BV{n) 2 


\\Tu-g\\l 2 ^n)+a J^\Du\, (1.4) 


variation, i.e., u € BV{Q) if and only if n G L^{Q) and 
\Du\ < oo. In the presence of impulse noise, e.g., 
salt-and-pepper noise or random-valued impulse noise, 
the above model usually does not yield a satisfactory 
restoration. In this context, a more successful approach, 
suggested in [HITlITSj, uses a non-smooth L^-data fi¬ 
delity term instead of the L^-data fidelity term in dEll), 
i.e., one considers 


min \\Tu- g\\Li(^f 2 )+a / \Du\, 
u€BV{Q) Jq 


(1.5) 


which we call the L^-TV model. In this paper, we are 
interested in both models, i.e., the L^-TV and the L^- 
TV model, and condense them into 


mm 

u£BV{n) 


^Jr{u;g) ■-'Hr{u-,g) + a J^\Du\^ ( 1 . 6 ) 


to obtain a combined model for removing Gaussian or 
impulsive noise, where T-Lriw^g) := 7 llT’n — for 

T — 1,2. Note, that instead of dm one can consider 
the equivalent problem 


min XHr{u:g)-\- / \Du\, 

EBVin) Jn 


(1.7) 


where A = ^ > 0. Other and different fidelity terms 
have been considered in connection with other type of 


risk estimator method (SURE) |571I551[TU] and its gen¬ 
eralizations the generalized cross-validation 

method (GGV) gSlIMlISTlITB] . the L-curve method gHl 
the discrepancy principle m, and the variational 
Bayes’ approach [^. Further parameter selection meth¬ 
ods for general inverse problems can be found for ex¬ 
ample in |^l42l[88ll8^ . 

Based on a training set of pairs {gk,Uk), for k = 
1 , 2,..., G N, where gk is the noisy observation and 
Uk represents the original image, for example in [TBIIMI 
EH bilevel optimization approaches have been presented 
to compute suitable scalar regularization parameters of 
the corresponding image model. Since in our setting we 
do not have a training set given, these approaches are 
not applicable here. 

Applying the discrepancy principle to estimate a in 
dm or A in dm, the image restoration problem can 
be formulated as a constrained optimization problem of 
the form 


min / \Du\ subject to (s.t.) g) = Br 

u^BV(n) J Q 

( 1 . 8 ) 

where Br '■= ■^|■!7| with Vr > 0 being here a constant 
depending on the underlying noise, r = 1,2, and |17| 
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(a) noisy image 


(b) over-smoothed reconstruction 


(c) over-fitted reconstruction 


Fig. 1.1 Reconstruction of an image corrupted by Gaussian white noise with a “relatively” large parameter a in (b) and a 
“relatively” small parameter oi in (c). 


denoting the volume of 17; see Section |2] for more de¬ 
tails. Note, that here we assume to know a-priori the 
noise level. In real applications this means that possibly 
in a first step a noise estimation has to be performed 
before the discrepancy principle may be used. However, 
in general it is easier to estimate the noise level than 
the regularization parameter |18| . 

The constrained minimization problem (11.81) is nat¬ 
urally linked to the unconstrained minimization prob¬ 
lem ini and accordingly to In particular, there 

exists a constant A > 0 such that the unconstrained 
problem dni) is equivalent to the constrained problem 
if T does not annihilate constant functions, i.e., 
T € £(L^(17)) is such that T • 1 = 1; see Section [2] 
for more details. Several methods based on the discrep¬ 
ancy principle and problem (HH) with T = 2 have been 
proposed in the literature, see for example [SUHHIIIM] 
and references therein, while not so much attention has 
been given to the case r = 1, see for example 


1.3 Spatially adaptive parameter 

Note, that a scalar regularization parameter might not 
be the best choice for every image restoration problem, 
since images usually have large homogeneous regions as 
well as parts with a lot of details. Actually it seems ob¬ 
vious that a should be small, or A should be large, in 
parts with small features in order to preserve the de¬ 
tails. On the contrary a should be large, or A should be 
small, in homogeneous parts to remove noise consider¬ 
able. With such a choice of a spatially varying weight we 
expect better reconstructions than with a globally con¬ 
stant parameter, as demonstrated for example in [3Zl 
[S5] . This motivated to consider multi-scale total varia¬ 
tion models with spatially varying parameters initially 
suggested in [SD] . The multi-scale version of (HU) reads 


as 

min 'Hr{u;g)+ / a{x)\Du\ (1-9) 

u^BV(n) Jq 

while for (El) one writes 

min — / \(x)\Tu — gVdx + [ \Du\, (1-10) 

ueBV{n) T Jq Jq 

and in the sequel we refer to da and (II.101) as the 
multi-scale L'^-TV model. 

In [53] the influence of the scale of an image fea¬ 
ture on the choice of a is studied and the obtained 
observations were later used in [55] to construct an up¬ 
dating scheme of a. Based on (11.101) in |5] a piecewise 
constant function A, where the pieces are defined by 
a partitioning of the image due to a pre-segmentation, 
is determined. In particular, for each segment a scalar 
Ai, i = 1,..., #pieces is computed by Uzawa’s method 

m- 

Later it was noticed that stable choices of A respec¬ 
tively a should incorporate statistical properties of the 
noise. In this vein, in [21 1571111] for the problem (11.101) 
automated update rules for A based on statistics of lo¬ 
cal constraints were proposed. In m a two level ap¬ 
proach for variational denoising is considered, where in 
the first level noise and relevant texture are isolated in 
order to compute local constraints based on local vari¬ 
ance estimation. In the second level a gradient descent 
method and an update formula for X(x) derived from 
the Euler-Lagrange equation is utilized. An adaptation 
of this approach to multiplicative noise can be found in 
[55] . For convolution type of problems in [2] based on 
an estimate of the noise variance for each pixel an au¬ 
tomatic updating scheme of A using Uzawa’s method is 
created. This approach is improved in [37] by determin¬ 
ing the fidelity weights due to the Gumbel statistic for 
the maximum of a finite number of random variables 
associated with localized image residuals and by incor¬ 
porating hierarchical image decompositions, proposed 
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in [HH1IR7] . to speed up the iterative parameter adjust¬ 
ment process. An adaptation of this approach to a total 
variation model with local constraints is studied in 
[58] . A different approach has been proposed in [85] for 
image denoising only, where non-local means [14] are 
used to create a non-local data fidelity term. While in 
all these approaches the adjustment of A relies on the 
output of T being a deteriorated image again, in [^ 
the method of m is adjusted to the situation where 
T is an orthogonal wavelet transform or Fourier trans¬ 
form. Very recently also bilevel optimisation approaches 
are considered for computing spatially adaptive weights 
[M1IM1I57] . 


1.4 Contribution 

Our first contribution of this paper is to present a meth¬ 
od which automatically computes the regularization pa¬ 
rameter a in (ira based on m for r = 1 as well as 
for T = 2. Our approach is motivated by the parameter 
selection algorithm presented in m, which was origi¬ 
nally introduced for L^-TV image denoising only, i.e., 
when T = I, where / denotes the identity operator. In 
this setting the algorithm in [18] is shown to converge 
to a parameter a* such that the corresponding mini- 
mizer Ua» of is also a solution of (11.81) . The proof 
relies on the non-increase of the function a i—>■ • 

However, this important property does not hold for op¬ 
erators T I in general. Nevertheless, we generalize 
the algorithm from [TB] to problems of the type (HD) 
for T = 1, 2 and for general linear bounded operators T, 
e.g., T might be a convolution type of operator. Utiliz¬ 
ing an appropriate update of a, which is different than 
the one used in [TB], we are able to show analytically 
and numerically that our approach indeed converges to 
the desired regularization parameter. Further, besides 
the general applicability of our proposed method it even 
possesses advantages for the case r = 2 and T = I over 
the algorithm from [18] with respect to convergence. 
More precisely, in our numerics it turned out that our 
proposed method always needs less or at least the same 
number of iterations as the algorithm from m till ter¬ 
mination. 

Motivated by multi-scale total variation minimiza¬ 
tion, the second contribution of this paper is concerned 
with the automated selection of a suitable spatially 
varying a for the optimization problem in (II.9|) for 
r = 1,2. Based on our considerations for an automatic 
scalar regularization parameter selection, we present al¬ 
gorithms where the adjustment of a locally varying a 
is fully automatic. Differently to the scalar case the ad¬ 
justment of a is now based on local constraints, sim¬ 


ilarly as already considered for example in dEZlEH]. 
However, our approach differs significantly from these 
previous works, where problem (I1.10|) is considered and 
Uzawa’s method or an Uzawa-like method is utilized for 
the update of the spatially varying parameter. Note, 
that in Uzawa’s method an additional parameter has 
to be introduced and chosen accordingly. We propose 
an update-scheme of a which does not need any addi¬ 
tional parameter and hence is not similar to Uzawa’s 
method. Moreover, differently to the approaches in m 
EH] where the initial regularization parameter A > 0 has 
to be set sufficiently small, in our approach any initial 
a > 0 is allowed. In this sense is our algorithm even 
more general than the ones presented in [^1^ . 

1.5 Outline of the paper 

The remaining of the paper is organized as follows: 
In Section [2] we revisit and discuss the connection be¬ 
tween the constrained minimization problem dni) and 
the unconstrained optimization problem dUD. Section 
[3] is devoted to the automated scalar parameter selec¬ 
tion. In particular, we present our proposed method 
and analyze its convergence behavior. Based on local 
constraints we describe in Section [3] our new locally 
adapted total variation algorithm in detail. Algorithms 
for performing total variation minimization for spatially 
varying a are presented in Section [5] where also their 
convergence properties are studied. To demonstrate the 
performance of the new algorithms we present in Sec¬ 
tion m numerical experiments for image denoising and 
image deblurring. Finally, in Section [7] conclusions are 
drawn. 


2 Constrained versus unconstrained 
minimization problem 

In this section we discuss the connection between the 
unconstrained minimization problem dUD and the con¬ 
strained optimization problem (ESI). For this purpose 
we introduce the following basic terminology. Let V be 
a locally convex space, V its topological dual, and (•, •) 
the bilinear canonical pairing over V x VA The domain 
of a functional J : V —K U {-boo} is defined as the set 

Dom(J') := (u G V : J{v) < oo}. 

A functional is called lower semicontinuous (l.s.c) if 
for every weakly convergent subsequence ^ v we 
have 

liminf > J"(D). 

y(n) 
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For a convex functional J : V —>■ K U {+oo}, we define 
the subdifferential oi at v € V, as the set valued 
function dj'{v) = 0 if J"(u) = oo, and otherwise as 

dj{v) = {u* e V : {v\u-v)+J{v) < J(u) Vu S V}. 

For any operator T we denote by T* its adjoint and 
by C{L‘^{n)) we denote the space of linear and con¬ 
tinuous operators from to Moreover, ga 

describes the average value of the function g G L^{fl) 
in fl defined by ga := g{x) dx. 

Theorem 2.1 A ssume that T G does not 

annihilate constant functions, i.e., TIq 0, where 
lfi{x) = 1 for X € f2. Then the problem 

min / \Du\ s.t. T-Lriu] g) < Br (2.1) 

u(^BV{n)Ja 

has a solution for r = 1,2. 

Proof For a proof we refer the reader to m and m- 

□ 

Moreover, we have the following statement. 

Proposition 2.1 Assume that T G C{LA{fi)) is such 
that T • 1 = 1 and < \\g — Then prob¬ 

lem (EH) is equivalent to the constrained minimization 
problem (HHl) for r = 1,2. 

Proof For r = 2 the statement is shown in [20]. We 
state the proof for r = 1 by noting it follows similar 
arguments as for r = 2. Let m be a solution of EH- 
Note, that there exists u G BV (f?) such that u = u -\- 
ga- We consider now the continuous function /(s) = 
||T(s'u + go) - g||Li(j 7 ) for s G [0,1]. Note that /(I) = 
\\Tu - g\\LpQ) < vi\n\ and /(O) = \\Tgn - g\\L^(^a) = 
Wg - goh 1 ( 17 ) ^ since T • 1 = 1, and hence 

there exists some s G [0,1] such that /(s) = i^i|l7|. 
Set u' = su which satisfies \\Tu' — gW^HQ) = and 

[ \Du \ = s [ \Du\ < liminf [ \Dun\, 

Jn Jo n^ao 

where {un)n is a minimizing sequence of EH. Hence 
u' is a solution of EH. □ 

Now we are able to argue the equivalence of the 
problems EH and EH- 

Theorem 2.2 Let T G £(L^(i7)) be such that T-1 = 1 
and Vr\L2\ < ||5 — go\\[^T(^Qy Then there exists a > 0 
such that the constrained minimization problem EH is 
equivalent to the unconstrained problem EH ; i.e., u is 
a solution of EH if and only if u solves (EH- 


Proof For r = 2 the proof can be found in m Prop. 
2.1]. By similar arguments one can show the statement 
for T = 1 , which we state here. 

Set TZ{u) := \Du\ and 

-hoc if ||m - glliqr?) > 

0 if ||m - fflliRr?) < 

Notice, that TZ and Q are convex l.s.c functions and 
problem EH is equivalent to min^ TZ{u) 0(Tu). We 
have Dom(7^) = BV{L2) fl L'^{L2) and Dom(t/) = {u G 
L^(17) : Q(u) < -l-oo}. Since g G TDom(7^), there ex¬ 
ists u G Dom(7^) with IjTu — (/||ii(i 7 ) < v\\fl\l2. As 
T G C{LA{fI)) is continuous, Q oT \s continuous at u. 
Hence, by [321 Prop. 5.6, p. 26] we obtain 

d{n + go T){u) = dn{u) + d{g o t){u) 

for all u. Further, g is continuous at Tu, and hence by 
m Prop. 5.7, p. 27] we have for all u, 

d{goT){u) = T*dg{Tu) 

where dg{u) = { 0 } if ||m— g||ii(r 2 ) < vi\L2\ arrddg{u) = 
{ad{\\u-g\\L^(^o)),a> 0 } if \\u-g\\L^f^o) = vi\n\. 

If M is a solution of EH and hence of EH, then 

0 G a(7^ + a o T){u) = dn{u) + T*dg{Tu). 

Since any solution of (11.81) satisfies \\Tu — gWh^^o) = 
zzi|l7|, this shows that there exists an a > 0 such that 

G€dn{u)+aT*d{\\Tu-g\\mo))- 

Hence for this a > 0, u is a minimizer of the problem 

in (11.61) . 

Conversely, a minimizer u of EH with the above a 
is obviously a solution of (11.81) with \\Tu — g\\L^(o) = 
zzi|l7|. This concludes the proof. □ 

Note, that ||TitQ, — gWr'^^o) is (only) convex with 
respect to Tua, and hence a minimizer of EH is in 
general not unique even in the simple case when T = /, 
i.e., for two minimizers it],, and it^ in general we have 
Tm], 7 ^ Tu^. On the contrary, \\Tua —g\W 2 (^Q) is strictly 
convex with respect to Tua, i.e., for two minimizers u\ 
and u\ of p.4l) we have Tu\^ = Tit^. Moreover, the 
function a i—>■ 'Hi{ua]g) is in general not continuous 
[35] . while a >->■ 'H 2 {ua',g) indeed is continuous [30] . 
where Ua is a respective minimizer of (EH- Hence we 
have the following further properties: 

Lemma 2.1 Let Ua be a minimizer of (EH then a i—> 
'Hriua^g) is non-decreasing for t = 1,2. Moreover, 
a^'H 2 {ua-,g) maps R+ onto [Q,\\g - go\\\ 2 ^^Q)]- 

Proof For a proof see [3(71133] . □ 
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Proposition 2.2 If Ua^ is a minimizer of 


3 Automated scalar parameter selection 


£iu,ai) := ||m-511^2(^2) 
for 1 = 1 , 2 , then we have 

ll^Cti ^Ct 2 11^2(22) < C ||5 



“ 5r3llL2(r2) ■ 


with C := min 


^2 

012 — 0(1 


02 —Ol 

n 

r 

O 2 +O 1 

5 

O 2 +O 1 

; 


Proof By [71 Lemma 10.2] we have 

Uci2 |Il2(22) — {£{'^012 j Q!i) £{Uai j CTl)) 

Qfl ' ' Qfl 

ll'^ai Uct2 |Il2(22) — {£{Uai j 0 : 2 ) £{Ua2 j 0 ^ 2 )) ■ 

02 ' 02 


In order to find a suitable regularization parameter 
o > 0 of the minimization problem (IlHI) we consider the 
corresponding constrained optimization problem (11.81) . 
Throughout the paper we assume that T does not anni¬ 
hilate constant function, which guarantees the existence 
of a minimizer of the considered optimization problems; 
see Sectional We recall, that in the constraint of (lEHl) 
the value Br is defined as Br = where 12 ,- G R 

is a statistical value depending on the underlying noise 
and possibly on the original image. 


3.1 Statistical characterization of the noise 


Summing up these two inequalities yields 


1 

ai 


1 

0:2 


Uc, - u. 


02 IIL 2 ( 22 ) 


< 


~ ~ 5lli2(22) - ||Uai - 5|li2(r2)) 


which implies 


\Uai ^Q2|Il2(22) 
Q !2 ~ oil 


< 


ai + a2 


(ll'^o!2 9\\\‘2i^n) ll^ai 5llL2(r2)) ■ 


( 2 . 2 ) 


By the non-decrease and boundedness of the function 
a 7^2 (uq; 5 ), see Lemma I^TTl it follows 


-*a 2 llL 2 ( 22 ) < 


Q !2 — Oil 


0:1 + Q ;2 

On the other hand inequality (12.21) implies 


II 5 5r2|lL2(22) ■ (2-3) 


Let us characterize the noise corrupting the image in 
more details by making similar considerations as in m 
Section 2]. Note, that at any point x € f2 the contami¬ 
nated image g{x) = J\f(Tu){x) is a stochastic observa¬ 
tion, which depends on the underlying noise. Two im¬ 
portant measures to characterize noise are the expected 
absolute value and the variance, which we denote by vi 
and 1^2 respectively. For images contaminated by Gaus¬ 
sian white noise with standard deviation cr, we typically 
set r = 2 and V 2 = If the image is instead corrupted 
by impulse noise, then we set r = 1 and we have to 
choose III properly. In particular, for salt-and-pepper 
noise ni G [min{ri, r 2 }, max{ri, r 2 }], while for random¬ 
valued impulse noise ni should be a value in the interval 
[|, |], where we used that for any point a; G 12 we have 
Tu{x) G [0,1]; cf. [SS]. Here vi seems to be fixed, while 
actually vi depends on the true (unknown) image u. In 
particular, for salt-and-pepper noise the expected ab¬ 
solute value is given by 


{ II "^02 - gWriiO) + II - 5'llL2(r2)) , 

where we used the binomial formula of — = (a -I- 

6 ) (a — h) for a, 6 G M and the triangle inequality. Using 
Lemma 12.11 yields 


OL 2 ~ oil 
oil -\- 012 


||uq,j '^a2\\L^{n) 2 


a2 — Q!i 


II 5 9o\\L2^Qy (2.4) 


Qfi -|- Qf 2 

The assertion follows then from (12.31) and (12.41) . □ 


Remark 2.1 Without loss of generality let q ;2 > ai in 
Proposition 12.21 then we easily check that 


C = 



2 0)2 — 0(1 
>. Q;2+Q;i 


if 02 > foi, 
if 02 = fai, 
otherwise. 


yi{u) := r 2 - {r 2 - ri)-^ JjTu){x) dx (3.1) 

and for random-valued impulse noise we have 
Mu) - Mu)ix) + dx. 

(3.2) 

However, instead of considering the constraint Hr (u; g) = 
Br {u) in (11.81) , which results in a quite nonlinear prob¬ 
lem, in our numerics we choose a reference image and 
compute an approximate value Br. Since our proposed 
algorithms are of iterative nature (see APS- and pAPS- 
algorithm below), it makes sense to choose the current 
approximation as the reference image, i.e., the reference 
image changes during the iterations. Note, that for salt- 
and-pepper noise with ri = r 2 the expected absolute 
value becomes independent of u and hence vi = ri. 
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In case of Gaussian noise Vt and Br are independent 
of u too. Nevertheless, in order to keep the paper con¬ 
cise, in the sequel instead of Vt and Br we often write 
Vr{u) and Br{u), where u represents a reference image 
approximating u, even if the values may actually be 
independent from the image. 


3.2 Automated parameter selection strategy 


while if Hriuoio', g) > Br(uo,^) then there has to exist 
an an+i < an such that T-Lriuan+Bg) ^ 

This holds true by setting in every iteration 


Q^n+1 • — 


BrjUgJ 

'^r{Ua„]g) 


P 


(3.3) 


together with an appropriate choice of p > 0. In partic¬ 
ular, there exists always a p > 0 such that this condition 
is satisfied. 


In order to determine a suitable regularization param¬ 
eter a in |18] an algorithm for solving the constrained 
minimization problem (11.81) for T = I and r = 2 is 
proposed, i.e., in the presence of Gaussian noise with 
zero mean and standard deviation a. This algorithm 
relies on the fact that a i—)■ 'H 2 {ua]g) is non-decreasing, 
which leads to the following iterative procedure. 


Chambolle’s parameter selection 

Choose ao > 0 and set n := 0. 

(CPS): 

1) Compute 


Ua„ G arg min ||u - 511^2(^2) -f 2a„ 

uGBV{S2) ^ ' 

Jn 

2) Update a„+i := ««• 

3) Stop or set n := n - 1 - 1 and return to step 1). 


For the minimization of the optimization problem in 
step 1 ) in [IH] a method based on the dual formulation 
of the total variation is used. However, we note that 
any other algorithm for total variation minimization 
might be used for solving this minimization problem. 
The CPS-algorithm generates a sequence {ua„)n such 
that for n —>■ oo, ||ua„ — g\\L^{n) cr-\/|I 2 | and Uq,„ 

converges to the unique solution of with T = I 

and r = 2 [TH]. The proof relies on the fact that the 

function a —>■ non-increasing. Note, that 

this property does not hold in general for operators T ^ 
I. 

3.2.1 The p-adaptive algorithm 

We generalize now the CPS-algorithm to optimization 
problems of the type (11.81) for r = 1,2 and for gen¬ 
eral operators T. In order to keep or obtain appropri¬ 
ate convergence properties, we need the following two 
conditions to be satisfied. Firstly, the function a i—!> 
TLriua'^g) has to be monotonic, which is the case due 
to Lemma 12.11 Secondly, in each iteration n the pa¬ 
rameter an has to be updated such that {'Hr{ua„', g))n 
is monotonic and bounded by {Br{uar))n- More pre¬ 
cisely, if TLriuao', g) < Br(uag) then there has to exist 
an an+l ^ 0.n such that TLriUan+i', g) < 


Proposition 3.1 Assume \\g — gnW'i^T^^Q'^ > and 

an+i is defined as in (1531) . 

(i) If an > 0 such that Hrinar', g) = Briua„), then for 
all p €M. we have that TLr {uar+i \g) ^ iua„+i ) • 
(a) If an > 0 such that 0 < TLriuan', g) < Br{ua^), then 
there exist p > 0 with "Hr(««„+!(ff) < Br{ua^^.^^). 

(Hi) If an > 0 such that ’Hriuan', g) > Br{ua„), then 
there exist p > 0 with 'Hr(««„+!(ff) > Briua„^i). 

Proof The assertion immediately follows by noting that 
for p = 0 we have an+i = an. □ 

Taking these considerations into account, a general¬ 
ization of the CPS-algorithm can be formulated as the 
following p-adaptive automated parameter selection al¬ 
gorithm: 

pAPS-algorithm: Choose oq > 0, p := po > 0, 

and set n := 0 . 

1) Compute e argmin„gBy(r 2 ) JV(u; p) 

2) Update a„+i := ( 

ll-TiuanTg) > 0 3 ,nd continue with step 3). 

Otherwise increase a„, e.g., an ■= 10a„, and go 

to step 1 ). 

3) Compute G argmin„gBy(r 3 ) Jr{u;g) 

4) a) a Hriugg-, g) < Br{UaJ 

(i) if < BriUg^^J gO tO Step 

5) 

(ii) if 'Hriug„^i; g) > Hr(ua„+i), decrease p, 
e.g., set p := p/ 2 , and go to step 2 ) 

b) A Hriugg-, g) > Br{ugg) 

(i) if > BriUg^^J gO tO Step 

5) 

(ii) if TLriug^^.j^; g) < Br{ug^^^), decrease p, 
e.g., set p := p/ 2 , and go to step 2 ) 

5) Stop or set n := n -I- 1 and return to step 2). 

Note, that due to the dependency of an+i on p a 
proper p cannot be explicitly computed, but only iter¬ 
atively, as in the pAPS-algorithm. 

The initial po > 0 can be chosen arbitrarily. How¬ 
ever, we suggest to choose it sufficiently large in order 
to keep the number of iterations small. In particular in 
our numerical experiments in Section [5] we set po = 32, 
which seems large enough to us. 









8 


Andreas Langer 


Proposition 3.2 The p APS-algorithm generates mono¬ 
tone sequences {an)n and (’Hr(Ma„; 3))„ such that 
(■Hr(itQn;5))„ is bounded. Moreover, if TLT{uao',g) > 
Br{uao) or Briua) < ^\\g - goWl for all a > 0, then 
io:n)n is also bounded. 

Proof If TLriuoLo', g) > Br{uao), then by induction and 
Lemma 12.11 one shows that 0 < a„+i < and 0 < 
< ’Hrina^; g) for all n £ N. Consequently 
{an)n and ("Hr( uq.^;5 ))rt are monotonically decreasing 
and bounded. 

If'Hr (Mao! 5) < Br{uaa), due to Lemma [2d1 we have 
that 0 < q;„ < a„+i and Hr(wa„;3) < Hr(Ma„+i;5) for 
all n £ N and hence (an)„ and (Hr (na„; 5))ra are mono¬ 
tonically increasing. Since there exists B* > 0 such that 
'B.Tiua„',g) < Br{ua^) < B* for all n £ N, see Section 
13.11 g))n is also bounded. If we additionally 

assume that Br{ua) < yllff —5f2||r for all a > 0 and we 
set B* := maxaBr{ua), then Theorem 12.21 ensures the 
existence of an a* > 0 such that 'HT{ua»',g) = Bf. By 
Lemma 12.II it follows that < a* for all n £ N, since 
Hr(uQ,„;g) < Br{ua,f) < Bf. Hence, (a„)„ is bounded, 
which finishes the proof. □ 


9Hr(ua„;5) + oind \Dua„\; See [2H1 Prop. 5.6 -I- Eq. 
(5.21), p.26]. Consequently there exist Va„ G 9 \B>Ua„ \ 
such that —anVa„ £ 9Hr(MQ„;<7) with lim„_>oo= 
Va. By [79l Thm. 24.4, p. 233] we obtain that —ava £ 
9Hr(ua; g) with Va G d \Dua \ and hence 0 £ dffriua, a 
for n —>■ 00 . 

If Hr(uaoJ 9 ) > Br, then as above we can show by 
induction that Q!„ > a„+i and Hr (tian; 5 ) > 'HT{ua„+i', g) 
Br. Thus we deduce that (Hr(ua„; 5 ))n and («„)„ are 
non-increasing and bounded. Note, that there exists 
an Of* > 0 with Hr('aa*;<?) = Hr such that for any 
a < Of*, Hr(««;<?) < Hr. Hence if Q!„ < a*, then 
T~LT{uan] 9 ) ^ Hr. This implies, that Hr(ua„;<?) = Hr 
and an+i = an- The rest of the proof is identical to 
above. □ 

Remark 3.1 The adaptive choice of the value p in the 
pAPS-algorithm is fundamental for proving convergence 
in Theorem EH In particular, the value p is chosen in 
dependency of a, i.e., actually p = p(a), such that in 

the case of a constant Br the function a 1 —>■ - 

is non-increasing; cf. Fig. l6.2l aL 


Since any monotone and bounded sequence converges 
to a finite limit, also (a„)„ converges to a finite value 
if one of the assumptions in Proposition 13.21 holds. For 
constant Br we are even able to argue the convergence 
of the pAPS-algorithm to a solution of the constrained 
minimization problem m- 

Theorem 3.1 Assume that Br{u) = Br is a constant 
independent of u and ||(/— > zzr|I?|. Then the 

pAPS-algorithm generates a sequence {oLn)n such that 
lim„^,oo a„ = d > 0 withTLr{ua] g) = lim^^oo 'Hr{ua„;g) 
= Br and Ua^ Ua € argmin„gx Jriu, a) forn —>■ 00 . 

Proof Let us start with assuming that TLriuao', g) < 
Br. By induction, we show that a„ < a„+i andHr(na„;ff) 
< Hr(wa„+i; 5 ) < Hr. In particular, if T-Lr{ua,,\ g) < Br 
then a„+i = «« > where p > 0 such 

that 'Hr{uar,+i'j g) < Br', cf. pAPS-algorithm. Then by 
Lemma im it follows that 

Br{Uar.',g) < 'Hr{Uar,+i',g) < Hr. 

Note, that there exists an a* >0 with Hrina* ',g) = Br, 
see Theorem l2.2l such that for any a> a*. Hr (mq,; g) > 
Br', cf. Lemma o If Cln> CX*, then Rrinan'^ g) > Hr. 
Hence TLr{ua„', 9 ) = Br and a„-i-i = Q;„. Thus we de¬ 
duce that the sequences (Hr(ua„; p))n, and (an)n are 
non-decreasing and bounded. Consequently, there ex¬ 
ists an a such that lim„_,.oo an = a with TLriua'Tg) = 

Br. Let a = limji_,.oo Hr {ua^ ; g), then H = Hr {ua', g) = 
Hr. By the optimality of Ua.,, we have that 0 £ dJr{ua„,oin 


3.2.2 The non p-adaptive case 

A special case of the pAPS-algorithm accrues when the 
value p is not adapted in each iteration but set fixed. 
For the case p = 1 (fixed) we obtain the following au¬ 
tomated parameter selection algorithm. 
APS-algorithm: Choose oq > 0 and set n := 0. 

1) Compute Ua„ £ argmin„gBy(r 2 ) JV(u, «„) 

2) Update a„+i := if Hr{uan\g) > 0 

and continue with step 3). Otherwise increase 
an, e.g., an '.= lOa^, and go to step I). 

3) Stop or set n := n -I- 1 and return to step 1). 

Even in this case, although under certain assump¬ 
tions, we can immediately argue the convergence of this 
algorithm. 

Theorem 3.2 For a > 0 let Ua be a minimizer of 
Jr{u,a). Assume that Br(u) = Br is a constant in¬ 
dependent of u, the function a 1 —>■ is non¬ 

increasing, and Up — guiWr > ^^rll^l- Then the APS- 
algorithm generates a sequence (a„)„ C K+ such that 
lim„^.ooQ;r!, = d > 0 , lim„^oo Hr (««„; p) = Hr and 
Uan converges to Ua £ argmin„gBy(j 7 ) J(u,a) forn —>■ 
00 . 

Proof We only consider the case when Hr{uao',g) < 
Br by noting that the case Hr (uao j p) > Hr can be 
shown analogous. By induction, we can show that a„ < 
an+i and Hr{uaF:9) If Hr(Mct„+i;p) < Hr. More pre- 
)cisely, if Hriua„',9) < Br then an+i = H 7 Fn~g)^n > 









Automated Parameter Selection for Total Variation Minimization in Image Restoration 


9 


an and by Lemma 12.11 it follows that Hr ; g) < 
'Hriuctn+i: g)- Moreover, by the assumption that a i—>■ 
jg non-increasing we obtain 'Hriuoin+i'-i g) ^ 
g) = Br- That is, 


[ 31 ]. The unique minimizer Ua,k+i is given by Ua.fc-i-i = 
(/ — Ps^K){z{ua,k))^ where K is the closure of the set 

{div^ : e e |^(a;)| < 1 Va; G 12 } 


'BriUan^g) < 'Hr(Wa„+i;g) < Br- 
The rest of the proof is analog to the one of Theorem 

O □ 


Nothing is known about the convergence of the APS- 
algorithm, if Br{-) indeed depends on u and Br{ua^) is 
used instead of a fixed constant. In particular, in our nu¬ 
merics for some examples, in particular for the applica¬ 
tion of removing random-valued impulse noise with r = 
0.05, we even observe that starting from a certain iter¬ 
ation the sequence (an)n oscillates between two states, 
see Fig. 6.7(c) This behavior can be attributed to the 
fact that, for example, if Hr{ua^) < Briua^), then it 
is not guaranteed that also Hr{ua^^i) < 
which is essential for the convergence. 

The second assumption in the previous theorem, 
i.e., the non-increase of the function a i—> 'H^{ua,g) 
can be slightly loosened, since for the convergence of 
the APS-algorithm it is enough to demand the non¬ 
increase starting from a certain iteration h > 0. That 
is, if there exists a region U C R’*' where a i—>■ 
is non-increasing and C 17, then the algorithm 

converges; see Fig. 16.21 Analytically, this can be eas¬ 
ily shown via Theorem 13.21 by just considering a^ as 
the initial value of the algorithm. If r = 2, similar to 
the CPS-algorithm, we are able to show the following 
monotonicity property. 


Proposition 3.3 // there, exists a constant c > 0 such 
that ||T*(Tm - g)\\L^a) = c\\Tu - g\\L^O) for all u G 

then the function a i—>■ is non-increasing, 

where Ua is a minimizer of J 2 {u, a). 


Proof We start by replacing the functional J 2 by a fam¬ 
ily of surrogate functionals denoted by S and defined 
for u,a G X as 

5 1 

S{u,a) := J 2 {u,a) + -||u - - ^\\T{u - a)\\l 2 ^i 2 ) 

= (5||zi-z(a)|||2(f2)-f 2a [ \Du\g,T) 

J Q 

where 5 > ||T|p, z{a) := a — jT*(Ta — g), and ip is 
a function independent of u. It can be shown that the 
iteration 


Ma,o G A, = argmin5(M,Ua,fe), fc > 0 (3.4) 

U 

generates a sequence (ua,k)k which converges weakly for 
fc —>■ 00 to a minimizer Ua of 772(u, a), see for example 


and Pk{u) := argmin„gif ||w— v\\i^ 2 (^Qy, see [TH]. Then 
for fc —>■ 00, let us define 


/(a) := \\P^K{ziu^))\\^,^^^ 


T*{Tuo. - g) 


L^(n) 


Since ||T*(Tm - g)\\L'^{n) = c\\Tu - it follows 

that /(a) = ^\j2'H2{ua\ g)- The assertion follows by 
applying m Lemma 4.1], which is extendable to in¬ 
finite dimensions, to / and by noting that the non¬ 
increase of a 1—^ implies the non-increase of a 1—>■ 

OL 


We remark that for convolution type of operators 
the assumption of Proposition 13.31 does not hold in 
general. However, there exist several operators T, rele¬ 
vant in image processing, with the property ||T*(Tm — 
ff)llL2(r2) = llTu-g||i,2(j7). Such operators include T = 
I for image denoising, T = lu for image inpainting, 
where 1 d denotes the characteristic function of the do¬ 
main D C f2, and T = S o A, where S' is a sub¬ 
sampling operator and A is an analysis operator of 
a Fourier or orthogonal wavelet transform. The lat¬ 
ter type of operator is used for reconstructing signals 
from partial Fourier data m or in wavelet inpainting 
pS] . respectively. For all such operators the function 

a I—>■ jg non-increasing and hence by setting 

P = ^ fixed in the pAPS-algorithm or changing the 
update of a in the APS-algorithm to 


^n+1•— 



n: 


where B2 is a fixed constant, chosen according to (dU, 
we obtain in these situations a convergent algorithm. 

We emphasize once more, that in general the non¬ 
increase of the function a 1—>■ is not guaran¬ 

teed. Nevertheless, there exists always a constant p > 0 
such that a 1—>■ ig indeed non-increasing. 

For example, p = 5 for operators T with the property 

\\T*{Tu-g)\\L 2 {n) = \\Tu- g\\L2^Qy, cf. Propsition [331 

In particular, one easily checks the following result. 

Proposition 3.4 Let 0 < a < / 3 , and Ua and Ufs min- 
imizers of Jr {■, O') and 77 t (-,/ 3 ), respectively, for t = 
1,2. Then < i'P-r i’u-a.g)) p < 

_In In OL _ 

\n'HT{uj3\g)—\n'HT{uo,\g) ' 
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4 Locally constrained TV problem 


In order to enhance image details, while preserving ho¬ 
mogeneous regions, we formulate, as in EZllSH], a locally 
constrained optimization problem. That is, instead of 
considering (ira we formulate 


min / \Du\ s.t. / w{x,y)\Tu — g\^{y)dy < Vr 
u^BV(n)jQ Jq 

(4.1) 

for almost every x € f2, where w is a normalized filter, 
i.e., w € X 17), and w > 0 on 17 x 17 with 


/ / w{x,y)dydx = l 

Jn J n 

and 

[ [ w{x,y)\(l){y)\^dydx>e\\(l)\\l.f^f2) 

Jo Jo 


(4.2) 


for all ^ G L'^(17) and for some e > 0 independent of (j); 

cf. 1 ^1^ . 


4.1 Local filtering 

In practice for w we may use the mean filter together 
with a windowing technique, see for example [571158) . In 
order to explain the main idea we continue in a discrete 
setting. Let 17^ be a discrete image domain contain¬ 
ing Ni X N 2 pixels, Ni,N 2 G N, and by |17^| = N 1 N 2 
we denote the size of the discrete image (number of 
pixels). We approximate functions u by discrete func¬ 
tions, denoted by The considered functions spaces 
are X = ^and Y = X x X. In what follows for 

all G X we use the following norms 

ll^''ll. := \\n%^(o^) = ( E 

\xeO’' j 

for 1 < r < -1-00. Moreover we denote by v!q the average 
value of G V, i.e, Uq ■= 'l2xeOf^ u^{x). The dis¬ 
crete gradient : X —^ F and the discrete divergence 

div^ : F —>■ X are defined in a standard-way by forward 
and backward differences such that div^ := — (V^)*; 
see for example [TH1I551I551IH5] . With the above nota¬ 
tions and definitions the discretization of the general 
function in (im is given by 

Jr{u^,a)--Hr{u^)+R^{u’^) (4.3) 

where = i|lrV -/||;, r G {1,2}, : X ^ 

X is a bounded linear operator, a G and 

i?a(u'‘) := E ot{x)\X'^u^{x)\i2 (4.4) 

xeo^ 


with \y \12 = \Jyl^y\ for every y = (j/i,j/2) G R^. In 
the sequel if a is a scalar or a = 1 in (031), we write 
instead of Ra or i?i just aR or i?, respectively, i.e., 

R{u^) = E |V'‘w'‘(x)|p 

xeO>^ 

is the discrete total variation of u in 17^, and we write 
Er instead of E^- to indicate that a is constant. Intro¬ 
ducing some step-size h, then for h —>■ 0 (i.e. the number 
of pixels X1X2 goes to infinity) one can show, similar 
as for the case a = 1, that Ra F-converges to a\Du\; 
see [T51IH5] . 

We turn now to the locally constrained minimiza¬ 
tion problem, which is given in the discrete setting as 

min R(u^) s.t. SJ Ju^) < — for all Xi G . 

(4.5) 

Here Vr is a fixed constant and 

SE{u>^):=^ E ^l(r"«")(^M)-/(^M)r 

xs,teii,j 

denotes the local residual at Xij G with being 
some suitable set of pixels around Xij of size M^j-, i.e., 
Mij = \Xij\. For example, in [571I551I55] for Xij the set 

with a symmetric extension at the boundary and with 
(jj being odd is used. That is, 17“^ is a set of pixels in 
a w-by-w window centered at Xij, i.e., Mij = for 
all f, j, such that 17“^ ^ 17^ for Xi^j sufficiently close to 
dfJ. Additionally we denote by a set of pixels in a 
window centered at Xij without any extension at the 
boundary, i.e., 

: maxjl - (i,i),< (s,t) 

< min|^^^y^,(Xi -i,X2 -j)| j. 

Hence C 17^ for all Xij G 17^. Before we analyze 
the difference between 17“^ and with respect to the 
constrained minimization problem 031), we note that, 
since does not annihilate constant functions, the 
existence of a solution of (H3|) is guaranteed; see m 
Theorem 2] [^ Theorem 2]. 

In the following we set Br := -^117^1. 

Proposition 4.1 (i) If is a solution of (14.51) with 
Xij = then Hr{u^) < Br- 
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(ii) If u is a solution of (US) with lij = lifj, then 
Hriu’^) <Br. 


in Proposition I4.ir iii. then a solution ui of the locally 
constrained minimization problem dni) satisfies 


Proof (i) Since is a solution of and is a 
set of pixels in a w-by-w window, we have 

> 7 E = hM)- 

Here we used that due to the sum over i,j each 
element (pixel) in appears at most times. 
More precisely, any pixel-coordinate in the set := 
{(z,j) : min{t-l, j-1, Afi-t,iV 2 -j} > occurs 
exactly w^-times, while any other pixel-coordinate 
appears strictly less than w^-times. This shows the 
hrst statement. 

(ii) For a minimizer of 63 we obtain 

Br>Y.Sl^iu^) 

= ^ E = Hr{u% 

which concludes the proof. □ 


'Hr{ui-,g)<—\12\ and [ \Dui\ > [ \Dus\ 

Jn Jo 

where Ug is a solution of 63- 

From Proposition 14.21 and Remark |4.II we conclude, 
since R{Ug) < R{u^) and J^\Dus\ < J^\Dui\, that 
Ug and Us are smoother than and ui, respectively. 
Hence the solution of the locally constrained minimiza¬ 
tion problem is expected to preserve details better than 
the minimizer of the globally constrained optimization 
problem. Since noise can be interpreted as fine details, 
which we actually want to eliminate, this could also 
mean, that noise is possibly left in the image. 

4.2 Locally adaptive total variation algorithm 

Whenever depends on u problem 63 results in a 
quite nonlinear problem. Instead of considering nonlin¬ 
ear constraints we choose as in Section[3]a reference im¬ 
age u and compute an approximate Vt = Vt{u). Note, 
that in our discrete setting for salt-and-pepper noise we 
have now 

:= r2 - (r2 - n)^ ^ {T^u'^){x) 
and for random-valued impulse noise we have 

^ E ^ (^{T>^u'^){xf - {T'^u^){x) + 0 . 


Note, that illij = fifj then by Proposition H.ll a mini¬ 
mizer of 63 also satisfies the constraint of the problem 

min R{u^) s.t. Hr{u^) < B^ (4.6) 

(discrete version of (12.11) ') but is in general of course not 
a solution of 63- 

Proposition 4.2 Let lij = Ug be a minimizer 

of und Uj^ be u Tn%n%m'izer of (63; R{Ug) < 

R{uf). 

Proof Assume that R{u^) > R{u^)- Since u’l is a so¬ 
lution of (14.61) it satisfies the constraint Hr{u^) < Br- 
By Proposition 14.11 we also have Hr{u^) < Br- Since 
R{Ug) > R{ui), Ug is not the solution of (14.61) which is 
a contradiction. Hence, Riu'f) < R{u^)- □ 


In our below proposed locally adaptive algorithms we 
choose as a reference image the current approximation 
(see LATV- and pLATV-algorithm below), as also done 
in the pAPS- and APS-algorithm above. Then we are 
seeking for a solution such that SJj{u'^) is close to 

T 

We note, that for large a > 0 the minimization of 
(63 yields an over-smoothed restoration and the 
residual contains details, i.e., we expect Hriu’f) > Br- 
Hence, if > -^ we suppose that this is due 

to image details contained in the local residual image. 
In this situation we intend to decrease a in the local 
regions Rij- In particular, we define, similar as in m 
155] . the local quantity ffj by 


n 





if Sl.iu’f,) > 

otherwise. 


Remark A.l Proposition 14.11 and its consequence are Note, that < 1 for all i, j and hence we set 
not special properties of the discrete setting. Let the 
Hlter w in 63 be such that the inequality in (113 
becomes an equality with e = l/|f2|, as it is the case 


a{xij) := 


M,- 


rs.teli, 


rts 


Ct{Xs,t)- 


(4.7) 
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On the other hand, for small a > 0 we get an under¬ 
smoothed image which still contains noise, i.e., we 
expect Hr{u^) < Br- Analogously, if we 

suppose that there is still noise left outside the residual 
image in Xij. Hence we intend to increase a in the local 
regions lij by defining 


ruj _ 


SUO 

— otherwise. 


and setting a as in (HTl) . Notice, that now > 

1. These considerations lead to the following locally 
adapted total variation algorithm. 


LATV-algorithm: Choose ao > 0, p := po > 0, 
and set n := 0. 

1) Compute ujj S argmin^^hg^f a„) 

2) (a) If ) > Br{u’^ ), then set 


:= maxi SI 




(b) If < Briu'lJ, then set 

fi,j ■ = max I min | Slj {u’lj, 
3) Update 


,e 


an+l(Xij) anixs,t)- 


4) Stop or set n := n -I- I and return to step I). 


Here and below £ > 0 is a small constant (e.g., in 
our experiments we choose e = 10“^^) to ensure that 
/“j > 0, since it may happen that Slj{u^^) = 0. 

If Hr{u^g) > Bt{u^^), we stop the algorithm as 
soon as the residual Hr{u^^) < Br{ul^^) for the first 
time and set the desired locally varying a* = a„. If 
HriUoa) — algorithm as soon 

as the residual > Br{u^^) for the first time 

and set the desired locally varying a* = a„_i, since 




,) 


The LATV-algorithm has the following monotonic¬ 
ity properties with respect to (an)n- 


Proposition 4.3 Assume Xi j = 17“^ and let e > 0 
be sujficiently small. If ao > 0 such that Hriu^g) < 
^r{Uag), then the LATV-algorithm generates a seguence 
(<ari)n such that 


^ ^ On-t-l ^ ^ ^ rXji {Xj^j ) ■ 

i,3 


Proof By the same argument as in the proof of Propo¬ 
sition O we obtain 


0^n-\-l{Xij) — 

ij i,j 

>E 

^,3 


( ^ 0^n{Xs,t) \ 

V ) 

( anix,j) \ 

\ rPcu^ ;■ 


Note that i/r(-) is bounded from below, see Section [XT] 
Consequently there exists an £ > 0 such that ^ > e 
for any Then, since Hriu^g ) if Br{u^Y have by 
the LATV-algorithm that 


/“ := max { min SfAu’f), 


VriutJ 


< 




and hence ,(«n)(u,i)- 


□ 


Proposition 4.4 Let Xi^j = and s > 0 be suffi¬ 
ciently small. 

(i) If uq > 0 such that Hriu^g) > Br{u^Y> then the 
LATV-algorithm generates a seguence (q;„)„ such 
that 

(ii) If ao >0 such that Hriu^Y — Briu^g), then the 
LATV-algorithm generates a seguence (an)n such 
that 

i,3 i,3 


Proof (i) By the same argument as in the proof of Propo- 
sition imi and since ffj := max1> 

- — we obtain 


E! ^n+l{Xi^j) — E! 

i,3 i,3 


'^r{Uan)^ OiniXs,t) 


= E 




— ^ ^ anixij). 
i,3 


(ii) Since Vt{-) is bounded from below, see Section IXTl 
there exists an £ > 0 such that ^ > e for any 
u^. Hence 


fz.i :=max min SfAu^), 


MKJ 


,e 



r 
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and by the same arguments as above we get 


Oin+i 


anjXij) 






In contrast to the pAPS-algorithm in the LATV- 
algorithm the power p > 0 is not changed during the 
iterations and should be chosen sufficiently small, e.g., 
we set p = ^ in our experiments. Note, that small p only 
allow small changes of a in each iteration. In this way 
the algorithm is able the generate a function a* such 
that Hr{u^,) is very close to . On the contrary, 

small p have the drawback that the number of itera¬ 
tions till termination are kept large. Since the parame¬ 
ter p has to be chosen manually, the LATV-algorithm, 
at least in the spirit, seems to be similar to Uzawa’s 
method, where also a parameter has to be chosen. The 
proper choice of such a parameter might be complicated 
and hence we are desiring for an algorithm where we do 
not have to tune parameters manually. Because of this 
and motivated by the pAPS-algorithm we propose the 
following p adaptive algorithm: 


pLATV-algorithm: Choose ao > 0, p := po > 0, 
and set n := 0. 

0) Compute Ua„ € argmin„hgjf 
1) (a) If Hriu^J > then set 



(b) If then set 



2) Update 





3) Compute G argmin^^g^ On+i) 

4) (a) ifi/.«)<i?.«) 

(i) if go to step 5) 

(ii) if i?r(Ua„+i) > -S^(Wa„+i). decrease p, 
e.g., set p = p/10, and go to step 2) 

(b) if > 5r«o) 

(i) if go to step 5) 

(ii) if I?r(Ua„+i) < decrease p, 

e.g., set p = p/10, and go to step 2) 

5) Stop or set n := n + 1 and return to step 1). 


In our numerical experiments this algorithm is ter¬ 
minated as soon as \Hr{Ua„) — Br{u'^^)\ < 10“® and 
< Br{u^^). Additionally we stop iterating 
when p is less than machine precision, since then any¬ 
way no progress is to expect. Due to the adaptive choice 
of p we obtain a monotonic behavior of the sequence 

Proposition 4.5 The sequence (q!„)„ generated by the 
pLATV-algorithm is for any point x G L2 monotone. In 
particular, it is monotonically decreasing for ao such 
that > Br{u^^), and monotonically increasing 

for ao such that < Br{u’fJ. 

Proof For Hr{u^^) > Briu^^) we can show by induc¬ 
tion that by the pLATV-algorithm f!fj > and 

hence 1 > ^ for all n. Then by the definition of 

an+i it follows 



^ C^n{Xi^jf 


By similar arguments we obtain for ao with Ht(u^q) < 
BriT^ao) O^n-t-l ^ ^n{Xi,j^ for all Xi^j G 12. D 

We are aware of the fact that using as a ref¬ 
erence image in the LATV- and pLATV-algorithm to 
compute Vr may commit errors. However, we recall that 
for Gaussian noise we set V 2 = and for salt-and- 
pepper noise with ri = r2 we have vi = ri. In these 
cases Vt does not depend on the original image and 
hence we do not commit any error by computing Vt. 
For random-valued impulse noise corrupted images the 
situation is different and vi indeed depends on the true 
image. In this situation errors may be committed when 
is used as a reference image for calculating Vr] see 
Figs. 14.11 and 14.21 Hence, in order to improve the pro¬ 
posed algorithm, for such cases for future research it 
might be of interest to find the optimal reference image 
to obtain a good approximation of the real value Vt- 

In contrast to the SA-TV algorithm presented in 
[SZllSH], where the initial regularization parameter has 
to be chosen sufficiently small, in the LATV-algorithm 
as well as in the pLATV-algorithm the initial value ao 
can be chosen arbitrarily positive. However, in the case 
Hr{Uag) > Br{u’f^) we cannot guarantee in general 
that the solution Ua obtained by the pLATV-algorithm 
fulfills Hriu’f) < Brivlff), not even if Br{ ) is con¬ 
stant, due to the stopping criterion with respect to the 
power p. On the contrary, if < Br{u^^), then 

the pLATV-algorithm generates a sequence such 

that < Br{u^^) for all n and hence also for the 
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Iteration Iteration 



Fig. 4.1 Progress of ) of the LATV-algorithm with P = ^ and ao = 10 ^ for removing random-valued impulse noise 

with r = 0.3 (left), r = 0.1 (middle), r = 0.05 (right) from the cameraman-image (cf. Fig. \6.lY b)}. 





Fig. 4.2 Progress of ) of the pLATV-algorithm with Po = ^ and ao = 10 ^ for removing random-valued impulse 

noise with r = 0.3 (left), r = 0.1 (middle), r = 0.05 (right) from the cameraman-image (cf. Fig. \6.lV h)). 


solution of the algorithm. As a consequence we would 
wish to choose Oq > 0 such that ^^r(Uao) — ^riu^g), 
which may be realized by the following simple auto¬ 
mated procedure: 

Algorithm 1: Input: oq > 0 (arbitrary); 

1) Compute G argmin„hgx Jriu^iCto)- 

2) If i?T(Uao) > Briuag) decrease ao by setting 
OiQ = Ca^ao, with Cag G (0,1), and continue with 
step 1), otherwise stop and return ag. 


and present solution methods, first for the case = I 
and then for general linear bounded operators T^. 

5.1.1 An algorithm for image denoising 

If = I, then (15.211 becomes an image denoising prob¬ 
lem, i.e., the minimization problem 

min ||m^ - g^|| 2 -I-2ii„(u''). (5.3) 


For solving this problem we use the algorithm of Cham- 
5 Total variation minimization bolle and Pock |22) . which leads to the following itera- 

In this section we are concerned with developing numer¬ 
ical methods for computing a minimizer of the discrete 
multi-scale L'^-TV model, i.e., 

min { Jt-(u^, a) := iJr(u^)-l-i?a(M^)}. (5.1) 


5.1 L^-TV minimization 

Here we consider the case r = 2, i.e., the minimization 
problem 

min i||TV-g'‘||i + i?„(w'^), (5.2) 

u'*ex z 


tive scheme: 


Chambolle-Pock algorithm: Initialize r, cr > 0, 

9 G [0,1], (pqjUq) gY X X, set Ug = Ug, and set 
n = 0. 

1 . Compute 

Pnix) + CrY'^Unix) 

Pn+lv^/ r 'I 5 

|^|Pn(^) + crV'*u()(a:)|, l| 

for all X G . 

2. Compute u„+i =- T+r -■ 

3. Set -f - <). 

4. Stop or set n := n -I- 1 and return to step 1). 
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In our numerical experiments we choose 9 = 1. 
In particular, in m it is shown that for 9 = 1 and 
rcrllV^lP < 1 the algorithm converges. 

5.1.2 An algorithm for linear bounded operators 

Assume, that is a linear bounded operator from X to 
X, different to the identity I. Then instead of minimiz¬ 
ing (15.21) directly, we introduce the surrogate functional 

5(u^ a'^) : = 

= - z{a’^)\\l + + i;{a\g\T'^), 

(5.4) 

with a’^,u^ G AT, z{a^) = - /), 

Ip a function independent of and where we assume 
6 > ||T^|p; see [?n[5^ . Note that 

min S{u^,a^) ^ min \\u^ - z{a^)\\l + 2R^{u^) 
'u'‘gx «ex ^ 

and hence to obtain a minimizer amounts to solve a 
minimization problem of the type o and can be 
solved as described in Section [5.1.11 Then an approxi¬ 
mate solution of dEll) can be computed by the follow¬ 
ing iterative algorithm: Choose Uq G X and iterate for 
n > 0 

= arg min S{u'^,u'P). (5.5) 

u''ex 

For scalar a it is shown in pSlI^TlfS^ that this iter¬ 
ative procedure generates a sequence {upfjn which con¬ 
verges to a minimizer of (Ol) . This convergence prop¬ 
erty can be easily extended to our non-scalar case yield¬ 
ing the following result. 

Theorem 5.1 For a : 17 —>■ R+ the scheme in (ESD 
generates a sequence {upP)n, which converges to a solu¬ 
tion of (El for any initial choice of Uq G X. 

Proof A proof can be accomplished analogue to [35]. 

□ 

5.2 An algorithm for L^-TV minimization 
The computation of a minimizer of 

min ||TV-/||i+ii„(u''), (5.6) 

is due to the non-smooth £^-term in general more com¬ 
plicated than obtaining a solution of the L^-TV model. 


Here we suggest to employ a trick, proposed in E for 
L^-TV minimization problems with a scalar regulariza¬ 
tion parameter, to solve (15.611 in two steps. In partic¬ 
ular, we substitute the argument of the £^-norm by a 
new variable v, penalize the functional by an L^-term, 
which should keep the difference between v and Tu — g 
small, and minimize with respect to v and u. That is, 
we replace the original minimization (15.61) by 

min - v^\\l + R^{u'^), (5.7) 

where 7 > 0 is small, so that we have g^ ~ — v^. 

Actually, it can be shown that (15.711 converges to (15.611 
as 7 —>■ 0. In our experiments we actually choose 7 = 
10“^. This leads to the following alternating algorithm. 


Lb 

-TVq, algorithm: Initialize a > 0, 

Uq G X and 

set 

n := 0. 



1) 

Compute 




= arg 

min \\v% + ^\\T>^u>P 

CN CN 

1 

1 

2) 

Compute 




G arg 

min ±\\T^u’^-g^-vf^^\\l+R^(u'^) 
u^^X 27 

3) 

Stop or set 

n := n -|- 1 and return to step 1). 


The minimizer in step 1) of the L^-TVq, algo¬ 
rithm can be easily computed via a soft-thresholding, 
i.e., = ST(T''uO - g\ 7), where 


(g'^{x)--i if/(a:)> 7, 

ST(g'',7)(a;) = < 0 if |g'‘(a;)| < 7, 

[g^{x)-\-'y iig'^{x)<-'y 

for all x G 17^. The minimization problem in step 2) is 
equivalent to 

arg min h\T’^u^ - g'^ - + R^^iu^) (5.8) 

Z 

and hence is of the type (ES). Thus an approximate 
solution of (ESI) can be computed as described above; 
see Section EH 

Theorem 5.2 The sequence {Un,v!f)n generated by the 
L^-TVa algorithm converges to a minimizer of EZD- 

Proof The statement can be shown analogue to E ■ Id 
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5.3 A primal-dual method for L^-TV minimization 

For solving (HU) with r = 1 we suggest, alternatively 
to the above method, to use the primal-dual method 
of [S5] adapted to our setting, where a Huber regular- 
isation of the gradient of u is considered; see [35] for 
more details. Denoting by ft a corresponding solution 
of the primal problem and p the solution of the asso¬ 
ciated dual problem, the optimality conditions due to 
the Fenchel theorem [39] are given by 

— divp(a;) = —kAu{x) + — - T*{Tu{x) — g{x)) 

P “f P- 

p Tu{x) - g{x) 

P + fi ma.x{/3,\Tu{x) - g{x)\} 

—p(a:) = — Vm(x) if |p(x)|;2 < a{x) 

7 

-p(x) = a{x) I if |p(x)|i2 = a(x), 

\Vu(x)\i 2 

for all X € f2, where k, /3, /i, and 7 are fixed positive 
constants. The latter two conditions can be summarized 
to -p(x) = ,nax{7°irih|vSDI,2} - Then setting q = -p 

and V = maxt/i^Ta-gl} to the following system of 

equation: 


by means of the peak signal-to-noise-ratio (PSNR) [llj . 
which is widely used as an image quality assessment 
measure, and the structural similarity measure (MSSIM) 
[90], which relates to perceived visual quality better 
than PSNR. When an approximate solution of the L^- 
TV model is calculated, we also compare the restora¬ 
tions by the mean absolute error (MAE), which is an 
L^-based measure defined as 

MAE = ||m - uWl^o), 

where u denotes the true image and u represents the ob¬ 
tained restoration. In general, when comparing PSNR 
and MSSIM, large values indicate better reconstruction 
than smaller values, while the smaller MAE becomes 
the better the reconstruction results are. 

Whenever an image is corrupted by Gaussian noise 
we compute a solution by means of the (multi-scale) L^- 
TV model, while for images containing impulsive noise 
the (multi-scale) L^-TV model is always considered. 

In our numerical experiments the CPS-, APS-, and 
pAPS-algorithm are terminated as soon as 

l^r (wq^ ; 5 ) — Hr (riQ,„) I _ 1 n-5 


0 = — max{/3, \Tu{x) — g{x)\}v + Tu{x) — g{x) 

0 = divq(x) -I- kAu{x) — — T*{Tu{x) — g{x)) 

P T /i 


0 = max{7a(x), |VM(x)|;2}q(x) — a{x)Vu{x) 

(5.9) 

for all X G 17. This system can be solved efficiently by 
a semi-smooth Newton algorithm; see Appendix [A] for 
a description of the method and for the choice of the 
parameters k, j3, /i, and 7. 

Note, that different algorithms presented in the lit¬ 
erature can also be adjusted to the case of a locally 
varying regularization parameter, such as [laiiiiEH]- 
However, it is not the scope of this paper to compare 
different algorithms in order to detect the most efficient 
one, although this is an interesting research topic in its 
own right. 


or the norm of the difference of two successive iterates 
an and an+i drops below the threshold Cq, = 10“^°, 
i.e., ||a„ — a„+i|| < Cq. The latter stopping criterion 
is used to terminate the algorithms if (q:„)„ stagnates 
and only very little progress is to expect. In fact, if our 
algorithm converges at least linearly, i.e., if there exists 
an Ea G (0,1) and an m > 0 such that for all n > to we 
have ||a„+i—CTooll < £c(||a„ —OfooH, the second stopping 
criterion at least ensures that the distance between our 
obtained result a and aoo is ||q! — ccooH < 

6.1 Automatic scalar parameter selection 

For automatically selecting the scalar parameter a in 
(II.6p we presented in Section [5] the APS- and pAPS- 
algorithm. Here we compare their performance for im¬ 
age denoising and image deblurring. 

6.1.1 Gaussian noise removal 


6 Numerical experiments 

In the following we present numerical experiments for 
studying the behavior of the proposed algorithms (i,e., 
APS-, pAPS-, LATV-, pLATV-algorithm) with respect 
to its image restoration capabilities and its stability 
concerning the choice of the initial value uq- The per¬ 
formance of these methods is compared quantitatively 


For recovering images corrupted by Gaussian noise with 
mean zero and standard deviation a we minimize the 
functional in (II.6|) by setting r = 2 and T = I. Then 
H2 = ^|17| is a constant independent of u. The auto¬ 
matic selection of a suitable regularization parameter 
a is here performed by the CPS-, APS-, and pAPS- 
algorithm, where the contained minimization problem 
is solved by the method presented in Section 15. 1.1 1 We 
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recall, that by m Theorem 4] and Theorem 13.11 it is 
ensured that the CPS- and the pAPS-algorithm gener¬ 
ate sequences (an)n which converge to a such that Ua 
solves (HH). In particular, in the pAPS-algorithm the 
value p is chosen in dependency of a, i.e., p = p(a), 

such that a —>■ - ig non-increasing, see Fig. 

I6.2f al. This property is fundamental for obtaining con¬ 
vergence of this algorithm; see Theorem 13.11 For the 
APS-algorithm such a monotonic behavior is not guar¬ 
anteed and hence we cannot ensure its convergence. 
Nevertheless, if the APS-algorithm generates a’s such 
that the function a —> non-increasing, then 

it indeed converges to the desired solution, see Theo- 
rem l3.2l Unfortunately, the non-increase of the function 
a —!> '^•^(’^^3) always, see Fig. I6.2f b'). 

For a performance-comparison here we consider the 
phantom-image of size 256 x 256 pixels, see Fig. I6.1f ab 
corrupted only by Gaussian white noise with cr = 0.3. 
In the pAPS-algorithm we set po = 32. The behavior 
of the sequence (an)n for different initial ao, i.e., ag € 
{1,10“^, 10“^} is depicted in Fig. 16.31 We observe that 
all three methods for arbitrary ag converge to the same 
regularization parameter a and hence generate results 
with the same PSNR and MSSIM (i.e., PSNR= 19.84 
and MSSIM= 0.7989). Hence, in these experiments, de¬ 
spite the lack of theoretical convergence, also the APS- 
algorithm seems to converge to the desired solutions. 
We observe the same behavior for different a as well. 

Looking at the number of iterations needed till ter¬ 
mination, we observe from Table 16.11 that the APS- 
algorithm always needs significantly less iterations than 
the CPS-algorithm till termination. This is attributed 
to the different updates of a. Recall, that for a fixed a„ 


in the CPS-algorithm we set := 

while in the APS-algorithm the update is performed 


as a 


APS _ 

n+1 


T'H2(Ua„-,g) 


an- Note, that 


TW2(«Q„;g) — 


,if- 


if 


i22\n\ 

'R2('Ua„;g) 

^ 21^1 ^ 
T'H2(uc<„;g) — 


V2\0\ 


W 2 I 7 

1 . 


- 

Hence 


> 1 and 




we obtain a^ 


< a! 


r'H2{uc„ 
< 


CPS 


^ '^n+1 ^ 


.APS 

^n+1 


if 


if 


'22\S^\ 


T'H2(Ua 
< 


> 1 and 


;s) — 

1. That is, 


> a! 


CPS 


> a, 


APS 


n+1 ^ '^n+1 

in the APS-algorithm a 


r'H2{uc„-,g) — ^ 

changes more significantly in each iteration than in the 
CPS-algorithm, which leads to a faster convergence with 
respect to the number of iterations. Nevertheless, ex¬ 
actly this behavior allows the function a -+ 


'H2(Ua-,g) 


to increase which is responsible that the convergence of 
the APS-algorithm is not guaranteed in general. How¬ 
ever, in our experiments we observed that the function 


'H2{ua-,g) 


only increases in the first iterations, but 


non-increases (actually even decreases) afterwards, see 
Fig. l6.2f bb This is actually enough to guarantee conver¬ 
gence, as discussed in Section [31 since we can consider 






(b) cameraman 



(c) Barbara 



(d) lena 


Fig. 6.1 Original images. 


the solution of the last step in which the desired mono¬ 
tonic behavior is not fulfilled as a “new” initial value. 
Since from this point on the non-increase holds, we get 
convergence of the algorithm. 

The pAPS-algorithm is designed to ensure the non¬ 
increase of the function a -+ choos- 

Q! 
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(a) CPS-algorithm 



(b) APS-algorithm 



(c) pPAS-algorithm 


Fig. 6.3 Denoising of the phantom-image corrupted with Gaussian white noise with a =■ 0.03. 


Table 6.1 Number of iterations needed for the reconstruction of the phantom-image corrupted by Gaussian white noise with 
different standard deviations cr. In the pAPS-algorithm we set po = 32. 




a = 0.3 


<7 = 0.1 


a = 0.05 


a = 0.01 

ao 

CPS 

APS 

pAPS 

CPS 

APS 

pAPS 

CPS 

APS 

pAPS 

CPS 

APS 

pAPS 

1 

55 

25 

8 

47 

21 

21 

46 

20 

20 

47 

17 

47 

10-1 

42 

18 

18 

38 

17 

6 

44 

20 

20 

46 

18 

46 

10-2 

43 

27 

43 

40 

20 

40 

40 

18 

40 

39 

17 

6 

10-1* 

43 

31 

43 

41 

22 

41 

41 

19 

41 

41 

20 

41 

10-^ 

43 

35 

43 

41 

23 

41 

41 

22 

41 

41 

19 

41 


ing p{a) in each iteration accordingly, which is done 
by the algorithm automatically. If p{a) = p = 1/2 
in each iteration, then the pAPS-algorithm becomes 
the CPS-algorithm, as it happens sometimes in prac¬ 
tice (indicated by the same number of iterations in Ta¬ 
ble |6J]). Since the CPS-algorithm converges [18], the 
pAPS-algorithm always yields p > \I2. In particular, 
we observe that if the starting value ao is larger than 
the requested regularization parameter a, less itera¬ 
tion till termination are needed than with the CPS- 
algorithm. On the contrary, if oq is smaller than the 
desired a, p = 1/2 is chosen by the algorithm to en¬ 
sure the monotonicity. The obtained result of the pAPS- 
algorithm is independent on the choice of po as visible 
from Fig. 16.41 In this plot we also specify the number 
of iterations needed till termination. On the optimal 
choice of po with respect to the number of iterations, 
we conclude from Fig. 16.41 that po = 32 seems to do a 
good job, although the optimal value may depend on 
the noise-level. 

Similar behaviors as described above are also ob¬ 
served for denoising other and real images as well. 

6.1.2 Image deblurring 

Now, we consider the situation when an image is cor¬ 
rupted by some additive Gaussian noise and addition¬ 
ally blurred. Then the operator T is chosen according to 


the blurring kernel, which we assume here to be known. 
For testing the APS- and pAPS-algorithm in this case 
we take the cameraman-image of Fig. IG.ll bl. which is 
of size 256 x 256 pixels, blur it by a Gaussian blurring 
kernel of size 5x5 pixels and standard deviation 10 and 
additionally add some Gaussian white noise with vari¬ 
ance cr^. The minimization problem in the APS- and 
pAPS-algorithm is solved approximately by the algo¬ 
rithm in (15.51) . In Fig. 16.51 the progress of for differ¬ 
ent cr’s, i.e., a G {0.3,0.1,0.05}, and different oo’s, i.e., 
ao G (1,10“^, 10“^} are presented. In these tests both 
algorithms converge to the same regularization param¬ 
eter and minimizer. From the figure we observe, that 
the pAPS-algorithm needs much less iterations than the 
APS-algorithm till termination. This behavior might be 
attributed to the choice of the power p in the pAPS- 
algorithm, since we observe in all our experiments that 
p > 1 till termination. 

6.1.S Impulsive noise removal 

It has been demonstrated that for removing impulsive 
noise in images one should minimize the L^-TV model 
rather than the L^-TV model. Then for calculating a 
suitable regularization parameter a in the T^-TV model 
we use the APS- and pAPS-algorithm, in which the 
minimization problems are solved approximately by the 
L^-TVa-algorithm. Here, we consider the cameraman- 
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(a) a — 0.3; pAPS-algorithm 



(d) a = 0.3; APS-algorithm 



-do = 1 

-ao = 10-^ 

■ ao = 10 “^ 


-ao = 1 
-ao = 10-1 

■ ao = 10“^ 


1 

0.9 
0.8 

0.7^ ‘- 

0.6 I 

0.5 

0.4 

0.3 

0.2 

0.1 ..V.._ 

0^ -^^^^^^ 

5 10 15 20 25 30 

Iteration 

(b) cr — 0.1; pAPS-algorithm 



(e) (j = 0.1; APS-algorithm 


0.6 

0.5 

0.4 

0.3 



10 20 30 40 50 

Iteration 

(c) cr = 0.05; pAPS-algorithm 



(f) cr = 0.05; APS-algorithm 


Fig. 6.5 Reconstruction of the cameraman-image corrupted by Gaussian blurring kernel of size 5x5 and standard deviation 
10. In the pAPS-algorithm we set po = 32. 


image corrupted by salt-and-pepper noise or random¬ 
valued impulse noise with different noise-levels, i.e., ri = 
r2 G {0.3,0.1,0.05} and r G (0.3,0.1,0.05} respectively. 
The obtained results for different cto’s are depicted in 
Fig. 16. Gl and Fig. 16.71 For the removal of salt-and-pepper 
noise we observe from Fig. 16.61 similar behaviors of the 
APS- and pAPS-algorithm as above for removing Gaus¬ 
sian noise. In particular, both algorithms converge to 
the same regularization parameter. However, in many 
cases the APS-algorithm needs significantly less iter¬ 
ations than the pAPS-algorithm. These behaviors are 
also observed in Fig. 16.71 for removing random-valued 
impulse noise as long as the APS-algorithm finds a so¬ 
lution. In fact, for r = 0.05 it actually does not converge 


but oscillates as depicted in Fig. 6.7(c) 


spectively, and compute automatically a spatially vary¬ 
ing A based on a local variance estimation. However, 
as pointed out in [37ll58) . they only perform efficiently 
when the initial A is chosen sufficiently small, as we will 
do in our numerics. On the contrary, for the LATV- and 
pLATV-algorithm any positive initial ag is sufficient. 

For the comparison we consider four different im¬ 
ages, shown in Fig. 16.11 which are all of size 256 x 256 
pixels. In all our experiments for the SA-TV-algorithm 
we use lij = f2fj , see and we set the window-size 
to 11 X 11 pixels in the case of Gaussian noise and to 
21 X 21 pixels in case of impulse noise. For the LATV- 
and pLATV-algorithm we use the window-size oj = 11, 
if not otherwise specified, and choose Po = 


6.2 Locally adaptive total variation minimization 

In this section various experiments are presented to 
evaluate the performance of the LATV- and pLATV- 
algorithm presented in Section 01 Their performance is 
compared with the proposed pAPS-algorithm as well 
as with the SA-TV-algorithm introduced in [37] for L^- 
TV minimization and in |58j for L^-TV minimization. 
We recall that the SA-TV methods perform an approxi¬ 
mate solution for the optimization problem in (11.101) . re- 


6.3 Gaussian noise removal 

6.3.1 Dependency on the initial regularization 
parameter 

We start this section by investigating the stability of 
the SA-TV-, LATV-, and pLATV-algorithm with re¬ 
spect to the initial regularization parameter, i.e., Aq 
for the SA-TV-algorithm and ao for the other algo¬ 
rithms, by denoising the cameraman-image corrupted 
by Gaussian white noise with standard deviation a = 
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(d) ri = r 2 = 0.3; APS-algorithm 


(e) ri = r 2 =0.1; APS-algorithm 


(f) ri = r 2 = 0.05; APS-algorithm 


6.6 Denoising of the cameraman-image corrupted by salt-and-pepper noise. In the pAPS-algorithm we set po = 32. 
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Fig. 6.7 Denoising of the cameraman-image corrupted by random-valued impulse noise. In the pAPS-algorithm we set po = 
32. 









































































































Automated Parameter Selection for Total Variation Minimization in Image Restoration 


21 



(a) 



(b) 


Fig. 6.2 Denoising of the phantom-image corrupted with 
Gaussian white noise with a = 0.03. (a) Plot of the fune- 
tion a —>■ fhg pAPS-algorithm with oq = 10“^. 

(b) Plot of the funetion a —>■ APS-algorithm 

with ao = 10“^. The desired monotone behavior is observed 
after the second iteration (red part of the eurve). 


0.1. In this context we also compare the difference of 
the pLATV-algorithm with and without using Algo¬ 
rithm 1 for computing automatically an initial param¬ 
eter, where we set Ca^ = The minimization prob¬ 
lems contained in the LATV- and pLATV-algorithm 
are solved as described in Section 15.1.11 For compar¬ 
ison reasons we define the values PSNRdig := max^o 
PSNR(ao)~tninc(o PSNR(q;o) and MSSIMdis := max^j, 
MSSIM(ao) — nrinao MSSIM(q!o) to measure the varia¬ 
tion of the considered quality measures. Here PSNR(ao) 
and MSSIM(ao) are the PSNR and MSSIM values of 
the reconstructions, which are obtained from the con¬ 
sidered algorithms when the initial regularization pa¬ 
rameter is set to oo- From Table [62] we observe that 
the pLATV-algorithm with and without Algorithm 1 
are more stable with respect to the initial regulariza¬ 
tion parameter than the LATV-algorithm and the SA- 
TV-algorithm. This stable performance of the pLATV- 
algorithm is reasoned by the adaptivity of the value p, 
which allows the algorithm to reach the desired resid¬ 
ual (at least very closely) for any oq- As expected, the 
pLATV-algorithm with Algorithm 1 is even more stable 
with respect to ao than the pLATV-algorithm alone, 
since, due to Algorithm 1, the difference of the actu¬ 
ally used initial parameters in the pLATV-algorithm is 
rather small leading to very similar results. Note, that 
if ao is sufficiently small, then the pLATV-algorithm 
with and without Algorithm 1 coincide, see Table 16.21 
for ao G {10“^, 10“^, 10“"^}. Actually in the rest of our 
experiments we choose ao always so small that Algo¬ 
rithm 1 returns the inputted ao- 


0.4 
0.35 
0.3 
0.25 
a 0.2 
0.15 
0.1 
0.05 
0 

Fig. 6.4 Regularization parameter a obtained by the pAPS- 
algorithm with different po for denoising the phantom-image 
corrupted with Gaussian white noise for different a. 


6.3.2 Dependency on the local window 

In Table 16.31 we report on the performance-tests of the 
pLATV-algorithm with respect to the chosen type of 
window, i.e., lij = fffj and lij = Dfj. We observe 
that independently which type of window is used the 
algorithm Hnds nearly the same reconstruction. This 
may be attributed to the fact that the windows in the 
interior are the same for both types of window. Nev¬ 
ertheless, the boundaries are treated differently, which 
leads to different theoretical results, but seems not to 
have significant influence on the practical behavior. A 
similar behavior is observed for the LATV-algorithm, 
as the LATV- and pLATV-algorithm return nearly the 
same reconstructions as observed below in Table O 
Since for both types of windows nearly the same re¬ 
sults are obtained, in the rest of our experiments we 
limit ourselves to always set lij = in the LATV- 
and pLATV-algorithm. 

Next, we test the pLATV-algorithm for different val¬ 
ues of the window-size varying from 3 to 15. Fig. 16.81 
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Table 6.2 PSNR and MSSIM of the reconstruction of the eameraman-image corrupted by Gaussian white noise with standard 
deviation a = 0.1 via the LATV- and pLATV-algorithm with different ao and via the SA-TV-algorithm with different Aq. In 
the LATV- and pLATV-algorithm we use Xij = with window-size 11 x 11 pixels in the interior and we set po = 



SA-TV 

LATV 

pLATV 

pLATV with Algorithm 1 

ao/Ao 

PSNR 

MSSIM 

PSNR 

MSSIM 

PSNR 

MSSIM 

PSNR 

MSSIM 

1 

27.82 

0.8155 

27.44 

0.8258 

27.37 

0.8260 

27.37 

0.8168 

10-1 

27.77 

0.8123 

27.59 

0.8211 

27.41 

0.8189 

27.38 

0.8166 

10-2 

27.71 

0.8107 

27.39 

0.8167 

27.37 

0.8167 

27.37 

0.8167 

10-3 

27.42 

0.8007 

27.40 

0.8167 

27.38 

0.8168 

27.38 

0.8168 

10-4 

27.56 

0.7792 

27.40 

0.8168 

27.38 

0.8168 

27.38 

0.8168 

PSNRdiff 

0.39646 

0.20257 

0.044473 


0.012704 

MSSIMdiff 

0.036322 

0.0091312 

0.0092963 

0.00019843 


Table 6.3 PSNR and MSSIM of the reconstruction of different images corrupted by Gaussian white noise with standard 
deviation a via the pLATV-algorithm with ao = 10“^. 


Image 

cr 

pLATV 

PSNR 

with X = Q 
MSSIM 

pLATV 

PSNR 

with X = fl 
MSSIM 

cameraman 

0.3 

22.47 

0.6807 

22.47 

0.6809 


0.1 

27.38 

0.8168 

27.37 

0.8165 


0.05 

30.91 

0.8875 

30.92 

0.8875 


0.01 

40.69 

0.9735 

40.68 

0.9735 

lena 

0.3 

22.31 

0.5947 

22.30 

0.5950 


0.1 

26.85 

0.7447 

26.87 

0.7448 


0.05 

30.15 

0.8301 

30.15 

0.8300 


0.01 

39.69 

0.9699 

39.68 

0.9699 


shows the PSNR and MSSIM of the restoration of the 
cameraman-image degraded by different types of noise 
(i.e., Gaussian noise with a = 0.3 or cr = 0.1, salt-and- 
pepper noise with ri = r2 = 0.3 or ri = r2 = 0.1, or 
random-valued impulse noise with r = 0.3 or r = 0.1), 
where the pLATV-algorithm with ao = 10“^ and po = 
1/2 is used. We observe that the PSNR and MSSIM are 
varying only slightly with respect to changing window- 
size. However, in the case of Gaussian noise elimination 
the PSNR and MSSIM increases very slightly with in¬ 
creasing window-size, while in the case of impulse noise 
contamination such a behavior cannot be observed. In 
Fig. 16.81 we also specify the number of iterations needed 
till termination of the algorithm. From this we observe 
that a larger window-size results in most experiments 
in more iterations. 

6.3.3 Homogeneous noise 

Now we test the algorithms for different images cor¬ 
rupted by Gaussian noise with zero mean and differ¬ 
ent standard deviations a, i.e., tr G {0.3, 0.1, 0.05, 0.01}. 
The initial regularization parameter ao is set to 10“^ in 
the pAPS-, LATV-, and pLATV-algorithm. In the SA- 
TV-algorithm we choose Aq = 10“^, which seems suffi¬ 
ciently small. From Table l^^ we observe that all consid¬ 
ered algorithms behave very similar. However, for a G 


(0.1, 0.05, 0.01} the SA-TV-algorithm most of the times 
performs best with respect to PSNR and MSSIM, while 
sometimes the LATV- and pLATV-algorithm have larger 
PSNR and MSSIM. That is, looking at these quality 
measures a locally varying regularization weight is pre¬ 
ferred to a scalar one, as long as a is sufficiently small. 
In Fig. I6.9l we present the reconstructions obtained via 
the considered algorithms and we observe that the LATV- 
and pLATV-algorithm generate visually the best re¬ 
sults, while the result of the SA-TV-algorithm seems in 
some parts over-smoothed. For example, the very left 
tower in the SA-TV-reconstruction is completely van¬ 
ished. This object is in the other restorations still visi¬ 
ble. For large standard deviations, i.e. tr = 0.3, we ob¬ 
serve from Table [R4l that the SA-TV method performs 
clearly worse than the other methods, while the pAPS- 
algorithm usually has larger PSNR and the LATV- and 
pLATV-algorithm have larger MSSIM. Hence, when¬ 
ever the noise-level is too large and details are consid¬ 
erably lost due to noise, the locally adaptive methods 
are not able to improve the restoration quality. 

6.3.f Non-homogeneous noise 

For this experiment we consider the cameraman-image 
degraded by Gaussian white noise with variance = 
0.0025 in the whole domain fl except a rather small 
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Table 6.4 PSNR- and MSSIM-values of the reconstruction of different images corrupted by Gaussian white noise with 
standard deviation a via pAPS-, LATV-, pLATV-algorithm with chq = 10“'^ and SA-TV-algorithm with Aq = 10“^. In the 
LATV- and pLATV-algorithm we useXij = Oij with window-size 11 x 11 pixels in the interior and we set po = 


Image 

(7 

pAPS ( 
PSNR 

scalar cn) 
MSSIM 

SA 

PSNR 

-TV 

MSSIM 

L7 

PSNR 

iTV 

MSSIM 

pL 

PSNR 

ATV 

MSSIM 

phantom 

0.3 

19.84 

0.7989 

19.83 

0.8319 

20.35 

0.8411 

20.31 

0.8432 


0.1 

28.97 

0.9644 

28.97 

0.9648 

29.50 

0.9680 

29.50 

0.9680 


0.05 

34.97 

0.9887 

33.77 

0.9867 

35.51 

0.9882 

35.51 

0.9882 


0.01 

48.88 

0.9994 

47.38 

0.9987 

49.46 

0.9993 

49.53 

0.9993 

cameraman 

0.3 

22.62 

0.6911 

22.03 

0.6806 

22.47 

0.6807 

22.47 

0.6807 


0.1 

27.31 

0.8109 

27.56 

0.7792 

27.40 

0.8168 

27.38 

0.8168 


0.05 

30.75 

0.8788 

31.60 

0.8929 

30.95 

0.8878 

30.91 

0.8875 


0.01 

40.51 

0.9731 

40.92 

0.9649 

40.73 

0.9737 

40.69 

0.9735 

barbara 

0.3 

21.22 

0.5022 

19.78 

0.4470 

21.05 

0.5032 

21.05 

0.5032 


0.1 

24.70 

0.7145 

25.53 

0.7292 

24.93 

0.7278 

24.93 

0.7278 


0.05 

28.22 

0.8514 

29.94 

0.8801 

28.49 

0.8584 

28.49 

0.8584 


0.01 

38.91 

0.9791 

40.56 

0.9809 

39.08 

0.9788 

39.08 

0.9788 

lena 

0.3 

22.42 

0.5930 

21.09 

0.5474 

22.33 

0.5951 

22.31 

0.5947 


0.1 

26.84 

0.7393 

27.31 

0.7528 

26.85 

0.7447 

26.85 

0.7447 


0.05 

30.06 

0.8261 

30.92 

0.8385 

30.16 

0.8307 

30.15 

0.8301 


0.01 

39.62 

0.9685 

39.81 

0.9660 

39.76 

0.9708 

39.69 

0.9699 


area (highlighted in red in Fig. 6.10(a) I, denoted by 17, 
where the variance is 6 times larger, i.e., = 0.015 in 

this part. Since the noise-level is in this application not 
homogeneous, the pLATV-algorithm presented in Sec¬ 
tion 0] has to be adjusted to this situation accordingly. 
This can be done by making Vt (here t = 2) locally 
dependent and we write Vr = Vr{u){x) to stress the 
dependency on the true image u and on the location 
a; G 17 in the image. In particular, for our experiment 
we set i/2 = 0.015 in 17, while i/2 = 0.0025 in 17\17. Since 
Vr is now varying, we also have to adjust the definition 
of Br and Br to 

Br{u) := / Vr{u){x) dx and Br{u^) := ^ Vr{u^){x), 


Gll'* 


respectively for the continuous and discrete setting. Mak¬ 
ing these adaptations allows us to apply the pLATV- 
algorithm as well as the pAPS-algorithm to the appli¬ 
cation of removing non-uniform noise. 

The reconstructions obtained by the pAPS-algorithm 
(with Po = 32 and oq = 10“^) and by the pLATV- 
algorithm (with Po = 1/2 and oq = 10 “^) are shown 


in Figs. 6.10(b) and 6.10(c) respectively. Due to the 


adaptive choice of a, see Fig. 6.10(d) where light col¬ 


ors indicate a large value, the pLATV-algorithm is able 
to remove all the noise considerably, while the pAPS- 
algorithm returns a restoration, which still retains noise 
in 17. 


6.4 Deblurring and Gaussian noise removal 

The performance of the algorithms for restoring im¬ 
ages corrupted by Gaussian blur with blurring kernel 


of size 5x5 pixels and standard deviation 10 and ad¬ 
ditive Gaussian noise with standard deviation a is re¬ 
ported in Table 16.51 Here we observe that the LATV- 
as well as the pLATV-algorithm outperform the SA- 
TV-algorithm for nearly any example. This observation 
is also clearly visible in Fig. 16.111 where the SA-TV- 
algorithm produces a still blurred output. The pAPS- 
algorithm generates very similar reconstructions as the 
LATV- and pLATV-algorithm, which is also reflected 
by similar PSNR and MSSIM. Similarly as before, the 
pAPS-algorithm performs best when cr = 0.3, while 
for smaller a the LATV-algorithm has always the best 
PSNR. 


6.5 Impulse noise removal 

Since it turns out that the LATV- and pLATV-algorithm 
produce nearly the same output, here, we compare only 
our pAPS- and pLATV-algorithm for L^-TV minimiza¬ 
tion, with the SA-TV method introduced in m. where 
a semi-smooth Newton method is used to generate an 
estimate of the minimizer of (ll.lOj) . For the sake of a fair 
comparison an approximate solution of the minimiza¬ 
tion problem in the pLATV-algorithm is solved by the 
semi-smooth Newton method described in Appendix lAl 
For the SA-TV method we use the parameters sug¬ 
gested in and hence lij = fffj. Moreover, we set 
Ao = 0.2 in our experiments which seems sufficiently 
small. In Table and Table O we report on the 
results obtained by the pAPS-, SA-TV-, and pLATV- 
algorithm for restoring images corrupted by salt-and- 
pepper noise or random-valued impulse noise, respec¬ 
tively. While the pAPS- and pLATV-algorithm produce 
quite similar restorations for both type of noises, the 
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Table 6.5 PSNR- and MSSIM-values of the reconstruction of different images corrupted by Gaussian blur (blurring kernel 
of size 5x5 pixels with standard deviation 10) and additive Gaussian noise with standard deviation cr via pAPS-, LATV-, 
pLATV-algorithm with oq = 10“^ and SA-TV-algorithm with Aq = 10“^. In the LATV- and pLATV-algorithm we use 
Xij = Oij with window-size 11 x 11 pixels in the interior and set po = 


Image 

<7 

pAPS ( 
PSNR 

scalar cn) 
MSSIM 

SA 

PSNR 

-TV 

MSSIM 

L7 

PSNR 

iTV 

MSSIM 

pL 

PSNR 

ATV 

MSSIM 

phantom 

0.3 

16.23 

0.6958 

15.65 

0.6632 

16.32 

0.6995 

16.31 

0.6997 


0.1 

17.86 

0.7775 

17.31 

0.7442 

17.98 

0.7914 

17.97 

0.7923 


0.05 

18.89 

0.7784 

19.20 

0.8193 

19.49 

0.8343 

19.39 

0.8279 

cameraman 

0.3 

21.04 

0.6410 

19.26 

0.5990 

20.86 

0.6272 

20.86 

0.6272 


0.1 

23.11 

0.7175 

22.64 

0.6957 

23.17 

0.7157 

23.15 

0.7156 


0.05 

24.14 

0.7562 

23.75 

0.7393 

24.22 

0.7573 

24.21 

0.7570 

barbara 

0.3 

20.58 

0.4556 

18.95 

0.4314 

20.42 

0.4517 

20.42 

0.4515 


0.1 

22.16 

0.5589 

22.09 

0.5687 

22.16 

0.5597 

22.16 

0.5597 


0.05 

22.87 

0.6245 

22.88 

0.6268 

22.92 

0.6273 

22.90 

0.6255 

lena 

0.3 

21.75 

0.5542 

20.10 

0.5278 

21.71 

0.5529 

21.69 

0.5528 


0.1 

24.44 

0.6496 

24.39 

0.6574 

24.50 

0.6514 

24.49 

0.6510 


0.05 

25.83 

0.7047 

25.81 

0.7091 

25.92 

0.7066 

25.91 

0.7062 


Table 6.6 PSNR- and MSSIM-values of the reconstruction of different images corrupted by salt-and-pepper noise withri = r 2 
via pAPS-, pLATV-algorithm with oto = 10“^ and SA-TV-algorithm with Aq = 0.2 and window-size 21 x 21. In the pLATV- 
algorithm we use lij = Oij with window-size 11 x 11 pixels in the interior and we set po = §■ 


Image 

ri = r2 

PSNR 

^PS (scalar 
MSSIM 

a) 

MAE 

PSNR 

SA-TV 

MSSIM 

MAE 

PSNR 

pLATV 

MSSIM 

MAE 

phantom 

0.3 

14.48 

0.7040 

0.0519 

15.28 

0.6540 

0.0605 

14.50 

0.7053 

0.0519 


0.1 

18.39 

0.8412 

0.0214 

19.57 

0.8703 

0.0196 

18.63 

0.8610 

0.0202 


0.05 

21.61 

0.9257 

0.0103 

22.81 

0.9362 

0.0103 

21.55 

0.9327 

0.0104 

cameraman 

0.3 

21.60 

0.7269 

0.0343 

21.34 

0.6871 

0.0390 

21.59 

0.7271 

0.0344 


0.1 

25.49 

0.8822 

0.0155 

25.80 

0.8774 

0.0157 

25.57 

0.8844 

0.0154 


0.05 

28.80 

0.9389 

0.0087 

28.50 

0.9251 

0.0095 

27.99 

0.9263 

0.0095 

barbara 

0.3 

21.56 

0.6242 

0.0486 

20.54 

0.5889 

0.0537 

21.51 

0.6253 

0.0488 


0.1 

25.49 

0.8729 

0.0211 

25.27 

0.8650 

0.0202 

25.81 

0.8759 

0.0208 


0.05 

28.46 

0.9338 

0.0118 

27.90 

0.9313 

0.0110 

28.34 

0.9325 

0.0121 

lena 

0.3 

23.29 

0.6807 

0.0360 

22.61 

0.6397 

0.0404 

23.32 

0.6811 

0.0359 


0.1 

27.60 

0.8508 

0.0151 

27.78 

0.8459 

0.0152 

27.99 

0.8530 

0.0148 


0.05 

29.45 

0.8946 

0.0096 

29.74 

0.8863 

0.0100 

29.53 

0.8931 

0.0097 


SA-TV algorithm seems to be outperformed in most 
examples. For example, in Fig. l6.l21 we observe that the 
pAPS- and pLATV-algorithm remove the noise consid¬ 
erable while the solution of the SA-TV method still con¬ 
tains noise. On the contrary for the removal of random¬ 
valued impulse noise in Fig. IB.l.'ll we see that all three 
methods produce similar restorations. 


7 Conclusion and extensions 

For L^-TV and L^-TV minimization including convo¬ 
lution type of problems automatic parameter selection 
algorithms for scalar and locally dependent weights a 
are presented. In particular, we introduce the APS- 
and pAPS-algorithm for automatically determining a 
suitable scalar regularization parameter. While for the 
APS-algorithm its convergence only under some assump¬ 
tions is shown, the pAPS-algorithm is guaranteed to 
converge always. Besides the general applicability of 
these two algorithms they also possess a fast numeri¬ 
cal convergence in practice. 


In order to treat homogeneous regions differently 
than fine features in images, which promises a better 
reconstruction, cf. Proposition 14.21 and Remark l4.11 al¬ 
gorithms for automatically computing locally adapted 
weights a are proposed. These methods are much more 
stable with respect to the initial oq than the state-of- 
the-art SA-TV method. Moreover, while in the SA-TV- 
algorithm the initial Aq > 0 has to be chosen sufficiently 
small, in our proposed methods any arbitrary oq > 0 is 
allowed. Hence the LATV- and pLATV-algorithm are 
much more flexible with respect to the initialization. 
By numerical experiments it is shown that the recon¬ 
structions obtained by the newly introduced algorithms 
are similar with respect to image quality measure to 
the restorations obtained by the SA-TV algorithm. In 
the case of Gaussian noise removal (including deblur¬ 
ring) for sufficiently small noise-levels reconstructions 
obtained by locally varying weights seem to be qualita¬ 
tively better than results with scalar parameters. On 
the contrary, for removing impulse noise a spatially 
varying a or A is in general not always improving the 
restoration quality. 
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Table 6.7 PSNR- and MSSIM-values of the reconstruction of different images corrupted by random-valued impulse noise 
via pAPS-, pLATV-algorithm with oq = 10“^ and SA-TV-algorithm with Aq = 0.2 and window-size 21 x 21. In the pLATV- 
algorithm we useXij = Oij with window-size 11 x 11 pixels in the interior and we set po = ^ in the pLATV-algorithm. 


Image 

r 

pf 

PSNR 

^PS (scalar 
MSSIM 

a) 

MAE 

PSNR 

SA-TV 

MSSIM 

MAE 

PSNR 

pLATV 

MSSIM 

MAE 

phantom 

0.3 

17.83 

0.8120 

0.0317 

18.68 

0.8012 

0.0319 

18.19 

0.8303 

0.0305 


0.1 

22.46 

0.9273 

0.0113 

23.83 

0.9278 

0.0100 

22.58 

0.9328 

0.0112 


0.05 

25.55 

0.9636 

0.0058 

26.56 

0.9642 

0.0054 

25.45 

0.9665 

0.0057 

cameraman 

0.3 

24.87 

0.8337 

0.0213 

23.48 

0.7583 

0.0237 

24.19 

0.7887 

0.0234 


0.1 

29.33 

0.9359 

0.0087 

27.72 

0.9087 

0.0089 

28.60 

0.9204 

0.0093 


0.05 

31.46 

0.9603 

0.0053 

30.53 

0.9478 

0.0052 

30.84 

0.9442 

0.0058 

barbara 

0.3 

24.24 

0.8040 

0.0301 

23.96 

0.7977 

0.0280 

24.24 

0.7992 

0.0302 


0.1 

29.20 

0.9355 

0.0118 

28.60 

0.9327 

0.0101 

28.91 

0.9305 

0.0120 


0.05 

31.95 

0.9650 

0.0065 

30.65 

0.9578 

0.0059 

31.85 

0.9640 

0.0066 

lena 

0.3 

26.80 

0.8124 

0.0205 

24.74 

0.7560 

0.0236 

26.63 

0.8082 

0.0208 


0.1 

30.34 

0.8965 

0.0092 

28.96 

0.8833 

0.0095 

30.07 

0.8918 

0.0095 


0.05 

31.36 

0.9189 

0.0062 

30.42 

0.9180 

0.0059 

31.08 

0.9159 

0.0063 


For computing a minimizer of the respective multi¬ 
scale total variation model we present first and second 
order methods and show their convergence to a respec¬ 
tive minimizer. 

Although the proposed parameter selection algo¬ 
rithms are constructed to estimate the parameter a in 
(fTTH) and (HU, they can be easily adjusted to find a 
good candidate A in (HU and (II.IOL respectively, as 
well. 

Note, that the proposed parameter selection algo¬ 
rithms are not restricted to total variation minimiza¬ 
tion, but may be extended to other type of regularizers 
as well by imposing respective assumptions that guar¬ 
antee a minimizer of the considered optimization prob¬ 
lem. In order to obtain similar (convergence) results 
as presented in Section [3] and Section H] the considered 
regularizer should be convex, lower semi-continuous and 
one-homogeneous. In particular for proving convergence 
results as in Section [3] (cf. Theorem Id.II and Theorem 
[3U an equivalence of the penalized and correspond¬ 
ing constrained minimization problem, as in Theorem 
12.21 is needed. An example of such a regularizer for 
which the presented algorithms are extendable is the 
total generalized variation |13j . 
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A Semi-smooth Newton method for solving 

dSU 

A semi-smooth Newton algorithm for solving (ISU can be 
derived similar as in 1581 by means of vector-valued variables. 
Therefore let g^ £ where 

N = N 1 N 2 , denote the discrete image intensity, the dual 
variable, the spatially dependent regularization parameter, 
and the observed data vector, respectively. Correspondingly 
we define £ R 2 iVxiV discrete gradient operator. 


g jgJVxJV ^jjg cjiggrete Laplace operator, T^ £ R^^^ 
as the discrete operator, and (T^)* as the transpose of T^. 
Here | • |, max{-, •}, and sign(-) are understood for vectors in 
a componentwise sense. Moreover, we use the function [| • |] : 

R^^ ->• R2^ with [|v'"|]i = [|u'^|]i+iv = -b 

for 1 < i < A. 

For solving 115.9ll in every step of our Newton method we 
need to solve 








(A.l) 


where 


Ai = [D’^icN) - D’^iv^)xA,^D’^isign{T^u^ - g^))]T\ 

S’l=Tu’^-g^-D\mfsyvt 

- g^) 

P + T 




ejv G R^ is the identity vector, D*^{v) is a diagonal matrix 
with the vector v in its diagonal, = max{/3, |T^u^—g^|}, 
= max{7a'", [|V'^u^|]}, 


= D^it/sJ with {tpji 
XA.,^ = D^(t-^y with {t-^^i 


M^{v) 




with 


if {m/sji = P, 
else; 

if (m-yji = 
else; 


V = {v^,Vyy £ 


Since the diagonal matrices D^irnpy and are in¬ 

vertible, we eliminate and 5q from nil). which leads to 
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Fig. 6.8 Restoration of the cameraman-image eorrupted by 
different types of noise via the pLATV-method with different 
window-sizes. 



(a) pAPS (PSNR: 27.31; MSSIM: 
0.8109) 



(b) SA-TV (PSNR: 27.56; MSSIM: 
0.7792) 



(c) LATV (PSNR: 27.40; MSSIM: 
0.8170) 



(d) pLATV (PSNR: 27.38; MSSIM: 
0.8171) 


the following resulting system 


Fig. 6.9 Reeonstruetion of the cameraman-image eorrupted 
by Gaussian white noise with cr = 0.1. 


HkSu = fk 
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(a) noisy observation 



(b) pAPS (PSNR: 30.01; MSSIM: 
0.8551) 



(c) pLATV (PSNR: 30.85; MSSIM: 
0.8859) 




(d) Locally varying a 


(a) pAPS (PSNR: 23.11; MSSIM: 
0.7175) 



(b) SA-TV (PSNR: 22.64; MSSIM: 
0.6957) 



(c) LATV (PSNR: 23.17; MSSIM: 
0.7160) 



(d) pLATV (PSNR: 23.15; MSSIM: 
0.7160) 


Fig. 6.10 Reconstruction of the cameraman-image cor¬ 
rupted by Gaussian white noise with = 0.0025 except in 
the in (a) highlighted area where = 0.015. 


Fig. 6.11 Reconstruction of the cameraman-image cor¬ 
rupted by Gaussian white noise with cr = 0.1 and Gaussian 
blur. 
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(a) noisy image 


(b) pAPS (PSNR: 24.24; MSSIM 
0.8040) 


(c) SA-TV (PSNR: 23.96; MSSIM 
0.7977) 


(d) pLATV (PSNR: 24.24; MSSIM: 
0.7992) 


(a) noisy image 


(b) pAPS (PSNR: 21.56; MSSIM: 
0.6242) 


(c) SA-TV (PSNR: 20.54; MSSIM: 
0.5889) 


(d) pLATV (PSNR: 21.51; MSSIM: 
0.6253) 


Fig. 6.12 Reconstruction of the barbara-image corrupted by Fig. 6.13 Reconstruction of the barbara-image corrupted by 
salt-and-pepper noise with ri = r 2 = 0.3. random-valued impulse noise with r = 0.3. 
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where 

Hk := j-i^h 

/3 + M tJ- + P 

h ■■= 

M + P 

In general and hence Hk is not symmetric. In m it is 
shown that the matrix at the solution {Uk,v^,q^) = 
{u, V, q) is positive definite whenever 

[\qi\]i<{c.% and (|n^|), < 1 (A.2) 

for i = 1,. . ., AI. 

In case these two inequalities are not satisfied we project 
q^ and onto their feasible set, i.e., ((g^ )i, (gjl )i+Ar) is 
set to (a'*)imax{(a'‘)i, [kfe |]i}“^((gfe)i, (gfe)i+jv) and (v^)i 
is replaced by max{l, (|nj|)i}(v^)i. Then the modified sys¬ 
tem matrix, denoted by H^ is positive definite; see 1361 . As 
pointed out in 1581 we may use + ek-D^(ejv) with k = 0, 
Sk > 0 and sj, 4' 0 ^.s fc ^ oo instead of /t > 0 to obtain a pos¬ 
itive definite matrix. Then our semi-smooth Newton solver 
may be written as in 1581 : 

Semi-smooth Newton method: Initialize (ug,qg) G 

jjiV ^ jg2Ar Q 

1. Determine the active sets X. 4 , 3 ,, £ and XA^i_ £ 

jU2Arx2iV 

2. If (IA.2II is not satisfied, then compute H '^; otherwise 
set := Hk- 

3. Solve H'^ Su = fk for S^. 

4. Compute 5q by using 5^- 

5. Update := and := -h 5,. 

6. Stop or set fc := fc -|- 1 and continue with step 1). 

This algorithm converges at a superlinear rate, which fol¬ 
lows from standard theory; see [mg]. 

In our experiments we always choose re = 0, /3 = 10“®, 
7 = 10“^, and /i = 10®. 
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